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Abstract 


With  increased  computational  power  and  progress  in  numerical  methods  over  the 
past  several  decades,  Computational  Fluid  Dynamics  (CFD)  is  now  used  routinely 
as  a  powerful  tool  in  the  design  of  aircraft.  Current  production  CFD  codes  used  in 
the  aerospace  industry  are  usually  second  order  accurate.  High-order  methods  have 
the  potential  to  achieve  higher  accuracy  at  less  cost  than  low-order  methods.  This 
potential  has  been  demonstrated  conclusively  for  smooth  problems  in  the  latest  In¬ 
ternational  Workshops  on  High-Order  Methods.  For  non-smooth  problems,  solution 
based  hp-adaptation  offers  the  best  promise.  Adjoint-based  adaptive  methods  have 
the  capability  to  dynamically  distribute  computing  resources  to  areas  most  impor¬ 
tant  for  predicting  engineering  parameters,  such  as  lift  and  drag  coefficients.  The 
primary  objective  of  the  present  study  is  to  develop  robust  and  efficient  high-order 
CFD  methods  and  tools  for  the  compressible  Navier-Stokes  equations  that  can  pro¬ 
vide  engineering  accuracy  for  real  world  industry  problems.  This  report  outlines 
progresses  in  the  following  areas. 

The  flux  reconstruction  (FR)  or  the  correction  procedure  via  reconstruction 
(CPR)  method  used  in  this  work  is  a  high-order  differential  formulation.  We  develop 
a  parallel  adjoint-based  adaptive  CPR  solver  which  can  work  with  any  element-based 
error  estimate  and  handle  arbitrary  discretization  order  for  mixed  elements.  First, 
a  dual-consistent  discrete  form  of  the  CPR  method  is  derived.  Then,  an  efficient 
and  accurate  adjoint-based  error  estimation  for  the  CPR  method  is  developed  and 
its  accuracy  and  effectiveness  are  verified  for  the  linear  and  non-linear  partial  dif¬ 
ferential  equations  (PDE).  The  current  method  has  been  applied  to  aerodynamic 
problems.  Numerical  tests  show  that  significant  savings  in  the  number  of  DOFs  can 
be  achieved  through  the  adjoint-based  adaptation. 

The  usage  of  large-eddy  simulation  (LES)  methods  for  the  computation  of  tur¬ 
bulent  flows  has  increased  substantially  in  recent  years.  By  resolving  large  energetic 
scales,  LES  has  the  potential  to  exhibit  good  performance  especially  for  vortex 
dominated  or  massively  separated  flows.  Due  to  the  disparate  length  scales  in  LES, 
high-order  methods  are  preferred  for  their  high  accuracy.  Recently,  models  of  the 
sub-grid  scale  (SGS)  stress  with  the  high-order  FR/CPR  method  have  been  eval¬ 
uated  on  3D  turbulent  flows  and  the  ID  Burgers’  equation.  Preliminary  studies 
show  that  implicit  LES  (ILES),  which  does  not  involve  any  SGS  model,  has  the 
most  efficient  and  accurate  results.  In  addition,  a  mathematical  analysis  of  the  scale 
similarity  is  performed,  revealing  that  the  ratio  of  the  resolved  stress  to  the  SGS 
stress  is  7s,  where  7  is  the  ratio  of  the  second  filter  width  to  the  first  filter  width, 
under  the  assumption  of  small  filter  width. 

For  high-order  methods  to  be  effective,  they  must  be  paired  with  liigh-order 
meshes  that  include  high-order  representations  of  curved  boundaries.  A  user-friendly, 
GUI-based  software  named  meshCurve  is  developed  to  convert  linear  unstructured 
meshes  to  curved  high-order  meshes  .  Using  the  state-of-the-art  algorithms,  the  soft¬ 
ware  reconstructs  the  curved  geometry  from  the  linear  surface  mesh,  while  retaining 
sharp  feature  curves  .  The  upgrade  process  is  automatic  and  does  not  require  a  CAD 
file.  meshCurve  may  be  used  by  CFD  practitioners  and  researchers  to  easily  produce 
quality  liigh-order  meshes  from  existing  low-order  meshes.  The  CGNS  standard  is 
used  as  the  file  format  for  input  and  output. 
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1  Introduction 


1.1  Background 

With  continuous  progresses  in  numerical  methods  and  computer  hardware  over  the  past 
several  decades,  Computational  Fluid  Dynamics  (CFD)  is  now  used  routinely  as  a  pow¬ 
erful  tool  in  the  design  of  aircraft.  Current  production  CFD  codes  used  in  the  aerospace 
industry  are  usually  second  order  accurate.  High-order  methods  have  the  potential  to 
achieve  greater  accuracy  at  less  cost.  This  potential  has  been  demonstrated  conclusively 
for  smooth  problems  in  the  recent  several  International  Workshops  on  High-Order  Meth¬ 
ods  [107].  A  variety  of  high  order  methods  have  been  developed.  Refer  to  several  books 
[49,  63,  106]  and  reviews  [29,  105]  for  the  state-of-the-art  and  recent  progresses  in  the 
development  of  such  methods.  The  primary  objective  of  the  present  study  is  to  develop 
robust  and  efficient  high-order  CFD  software  for  compressible  Navier-Stokes  equations 
that  can  provide  engineering  accuracy  for  real  world  industry  problems. 

The  numerical  method  we  consider  is  a  nodal  differential  high-order  formulation  named 
the  correction  procedure  via  reconstruction  (CPR)  method  [54]  [108].  This  formulation  has 
some  remarkable  properties.  The  framework  is  easy  to  understand,  efficient  to  implement 
and  can  recover  several  well  known  methods  such  as  the  discontinuous  Galerkin  (DG) 
[85,  8,  9,  23,  24,  83,  111],  the  spectral  volume  method  (SV)  [104,  110,  72]  and  the  spectral 
difference  methods  (SD)[71,  64,  78,  97,  70].  For  recent  development  with  CPR,  interested 
readers  can  refer  to  [55,  109,  58,  17,  35,  36,  18,  56,  118]. 

For  non-smooth  problems,  solution  based  hp-adaptations  offer  the  best  promise.  Adap¬ 
tive  methods  have  the  capability  to  dynamically  distribute  computing  resources  to  de¬ 
sired  areas  of  the  computational  domain,  achieving  required  accuracy  at  minimal  cost 
[74,  20,  27,  53].  Because  they  can  ensure  reliability  and  increase  the  robustness  of  high- 
order  methods,  adaptive  methods  have  received  considerable  attention  in  the  high-order 
CFD  community  [48,  98,  99,  31,  102,  116,  30]. 

The  effectiveness  of  adaptive  methods  highly  depends  on  the  accuracy  of  the  error 
estimation.  There  are  at  least  three  major  types  of  adaptation  criteria:  gradient  or 
feature  based  Berger  and  Colella  [15],  Warren  et  al.  [112],  Barth  [7],  Harris  and  Wang 
[46],  residual-based  Ainsworth  and  Oden  [1],  Baker  [5],  Johnson  [61],  Shih  and  Qin  [93], 
Gao  and  Wang  [34],  Cagnone  and  Nadarajah  [17],  and  adjoint-based  Giles  and  Pierce 
[41],  Venditti  and  Darmofal  [99],  Hartmann  and  Houston  [48],  Venditti  and  Darmofal 
[98],  Becker  and  Rannacher  [13],  Giles  and  Pierce  [42],  Park  [82],  Venditti  and  Darmofal 
[100],  Fidkowski  and  Darmofal  [31],  Wang  and  Mavriplis  [101],  Leicht  and  Hartmann  [65], 
Li  et  al.  [67],  Yano  and  Darmofal  [117],  Ceze  and  Fidkowski  [21],  Heuristic  feature-based 
criteria  perform  refinements  around  some  unique  flow  features,  such  as  large  gradients  or 
strong  vorticity,  but  because  they  do  not  directly  relate  to  the  output  variables  of  interest, 
they  cannot  provide  universal  and  robust  error  estimation  Zhang  et  al.  [119],  Venditti  and 
Darmofal  [98].  The  residual-based  error  criteria  targets  the  elements  which  have  large 
discretization  error,  flagging  them  for  refinement.  The  locally  defined  element-wise  error 
may  lead  to  false  refinements  in  convection-dominated  problems.  The  dual-weighted 
residual  method  proposed  by  Becker  and  Rannacher  [12]  relates  a  specific  functional 
output  directly  to  the  local  residual  by  solving  an  additional  adjoint  equation.  It  can 
capture  the  error  propagation  effects  inherent  in  the  hyperbolic  equations.  This  kind  of 
adjoint-based  error  indicator  has  been  shown  very  effective  in  driving  an  hp-adaptation 
procedure  to  obtain  a  very  accurate  prediction  of  the  functional  outputs  [41,  42,  13,  16, 
117,  21,  102].  Recently,  Fidkowski  and  P.L  Roe  developed  a  new  error  indicator  based 
on  the  entropy  variables,  which  can  be  interpreted  as  the  dual  solution  of  the  output 
of  entropy  balance  on  the  whole  domain.  It  can  be  obtained  directly  from  the  state 
variables  without  solving  extra  adjoint  equations  and  has  been  successfully  applied  to 
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inviscid,  viscous  and  turbulence  flows  Ficlkowski  and  Roe  [32,  33]. 

The  adjoint  solution  is  required  for  the  error  estimate  and  output-based  adaptation. 
There  are  two  approaches  to  obtain  the  adjoint:  the  continuous  adjoint  and  the  discrete 
adjoint.  It  has  been  shown  that  the  discrete  adjoint  solution  leads  to  a  more  accurate  error 
estimation  for  the  fine  grid  functional,  while  the  continuous  adjoint  gives  a  better  output 
estimation  when  the  primal  and  adjoint  solutions  are  well  resolved  [68].  However,  the 
discrete  adjoint  solution  should  be  consistent  with  the  exact  adjoint  from  the  continuous 
adjoint  equation.  It  is  well  known  that  dual  consistency  can  significantly  impact  the  con¬ 
vergence  rate  of  both  the  primal  and  adjoint  approximations.  There  are  several  possible 
sources  of  dual  inconsistency  in  a  high-order  discretization.  A  dual-consistent  discretiza¬ 
tion  with  variational  forms,  such  as  the  finite  element  and  DG  methods,  were  examined 
for  the  Euler  and  Navier-Stokes  equations  in  the  literature  [41,  42,  47,  62,  90].  More 
recently,  adjoint-based  error  estimation  for  summation-by-parts  finite-difference  methods 
have  been  studied  [50,  51,  14].  However,  the  analysis  of  the  dual  consistency  for  compact 
high-order  differential-type  methods  appears  lacking.  This  is  one  of  the  focuses  of  the 
present  study. 

The  high-order  CPR  method  can  handle  arbitrary  solution  order  on  mixed  grids.  The 
marked  candidate  elements  for  adaptation  can  be  modified  by  enriching  its  solution  order 
or  subdividing  its  element  or  resizing  its  grid.  Thus,  the  ways  to  increase  the  discretization 
resolution  can  be  generally  classified  into  3  categories:  h-refinement,  r-refinement  and  p- 
refinement.  For  h-refinement,  subdivision  is  performed  locally  for  each  candidate  element 
to  increase  the  total  DOFs.  R-refinement  or  the  moving  mesh  method  keeps  the  total 
number  of  nodes  the  same  but  moves  the  location  of  the  grid  locally  or  globally  [53] .  With 
p-refinement,  the  local  degree  of  approximation  polynomial  is  modified.  The  moving  mesh 
method  with  curved  elements  in  3D  is  still  an  on-going  research  area.  In  this  work,  we 
only  modify  the  solution  polynomial  order  locally  or  subdivide  elements  hierarchically 
for  adaptations.  Intuitively,  h-refinement  should  be  applied  to  the  discontinuities  and 
p-enrichment  is  appropriate  in  smooth  flow  regions.  However,  the  optimal  choice  between 
h-  or  p-refinement  is  not  a  trivial  problem,  which  is  studied  in  the  previous  research 
[16,  52,  96,  37,  11,  21,  25,  63,  102], 

Reynolds-averaged  Navier-Stokes  (RANS)  methods  have  been  used  almost  exclusively 
for  the  computational  fluid  dynamics  (CFD)  analysis  of  practical  engineering  turbulent 
flows  for  the  past  several  decades.  In  RANS,  all  turbulent  fluid  dynamic  effects  are  re¬ 
placed  by  a  turbulence  model.  RANS-based  techniques  are  successfully  used  in  the  indus¬ 
try  for  many  cases.  However,  the  behavior  in  the  vortex  dominated  or  massively  separated 
flows  are  far  from  being  satisfactory.  On  the  other  end,  direct  numerical  simulation  (DNS) 
methods  resolve  all  turbulent  motions.  Without  the  influence  of  the  turbulence  modeling, 
it  gives  the  whole  spectrum  of  the  turbulent  flow.  But  DNS  will  still  remain  impractical 
for  its  high  cost  of  computational  power.  Large-eddy  simulation  (LES)  is  a  compromise 
of  these  two  ends.  In  LES,  large  energetic  scale  motions  are  resolved  while  the  small 
scale  motions  are  taken  care  of  by  the  SGS  models.  With  the  resolving  of  the  important 
scales,  the  solution  given  by  LES  is  expected  to  be  more  accurate  than  RANS,  but,  still 
affordable  by  the  computation  power  today. 

In  LES,  the  large  scales  and  small  scales  are  separated  by  a  low-pass  filter.  While 
the  large  scales  are  fully  resolved,  the  small  scales  are  believed  to  be  more  universal,  and 
thus  easier  to  model  than  the  large-scale  ones.  Many  SGS  models  have  been  developed  in 
the  last  four  decades.  We  focuses  on  five  of  them:  the  static  Smagorinsky  model  (SS)J. 
[57]  DK  [26],  the  dynamic  Smagorinsky  model  (DS)Germano  M  and  WH  [38],  the  scale- 
similarity  model  (SSM)Bardina  J  [6],  the  mixed  model  (MM)  Bardina  J  [6]  and  the  linear 
unified  RANS-LES  model  (LUM)  Harish  Gopalan  [45].  Among  explicit  models,  the  SS  is 
a  popular  one  because  of  its  simplicity.  The  effect  of  the  SGS  stress  upon  the  resolved 
scales  is  modeled  as  an  eddy  viscosity.  The  eddy  viscosity  is  expressed  in  the  mixing 
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length  form  with  a  dimensionless  empirical  coefficient.  However,  it  has  been  found  that 
the  empirical  coefficient  depends  on  the  flow.  It  also  adds  too  much  dissipation  to  the 
large  scale  motions  if  we  keep  the  coefficient  the  same  as  we  approach  wall  boundaries. 
To  resolve  these  deficiencies,  the  DS  model  was  developed  in  Germano  M  and  WH  [38]. 
In  the  DS  model,  the  coefficient  is  calculated  based  on  the  Germano  identity,  which 
involves  two  levels  of  filtering  and  relates  the  SGS  stress  to  the  resolved  stress.  The 
coefficient  is  locally  decided  and  no  longer  a  prescribed  constant,  and  it  goes  to  zero  as  a 
wall  boundary  is  approached.  The  DS  model  has  been  applied  to  a  large  variety  of  flow 
simulations  Piomelli  U  [84]Ghosal  S  [40]Akselvoll  K  [3]Wu  X  [113].  An  alternative  way 
to  model  the  SGS  stress  is  offered  by  the  SSM  Bardina  J  [6] .  As  the  name  indicates, 
it  assumes  similarity  between  two  scales  of  stresses,  the  resolved  stress  and  the  SGS 
stress.  Numerical  tests  showed  that  energy  accumulated  at  small  scales  with  this  model 
Meneveau  C  [80].  To  remedy  the  problem,  the  DS  was  added  to  dissipate  the  energy, 
which  led  to  the  MM.  Recently,  hybrid  RANS-LES  models  have  drawn  much  research 
interest.  They  combine  RANS  with  LES  so  that  in  the  near  wall  region,  the  RANS  model 
is  used,  while  LES  is  employed  in  the  outer  region.  These  models  have  demonstrated  good 
accuracy  with  reasonable  cost  when  compared  to  a  pure  LES  approach.  The  linear  unified 
RANS-LES  model  (LUM)  was  developed  in  Harish  Gopalan  [45].  In  the  present  work, 
we  also  evaluate  the  LUM  model.  Finally,  we  also  consider  the  monotone  integrated  LES 
Durbin  P.  A  [28]  or  implicit  LES  (ILES)  Grinstein  F.  F.  [44],  in  which  no  explicit  SGS 
model  is  used.  In  ILES,  the  numerical  algorithm  has  its  numerical  dissipation  which  serves 
as  the  SGS.  The  obvious  advantage  of  ILES  is  its  lower  computational  cost  compared  with 
the  conventional  SGS  models. 

As  high-order  methods  become  increasingly  mature  for  real  world  applications,  high- 
order  boundary  treatment  becomes  a  critical  issue.  In  mainstream  finite  volume  solvers, 
the  wall  boundary  of  the  geometry  is  usually  represented  by  piecewise  straight  segments 
or  planar  facets.  However,  this  linear  boundary  representation  is  far  from  sufficient  for 
high-order  methods.  When  the  error  from  the  geometry  dominate  the  discretization  error, 
the  benefit  of  high-order  methods  is  largely  lost.  For  high-order  methods  to  be  effective, 
they  must  be  paired  with  high-order  meshes  that  include  high-order  representations  of 
curved  boundaries  [77,  115]. 

Unfortunately,  algorithms  for  generating  high-order  meshes  are  in  short  supply,  and 
software  support  is  limited.  The  lack  of  adequate  high-order  mesh  generators  has  prompted 
researchers  at  the  1st  International  Workshop  of  High-Order  CFD  to  prioritize  high-order 
meshing  as  a  pacing  item  for  future  research  [103].  A  second,  related  issue,  is  the  need 
for  high-order  versions  of  existing  low-order  meshes.  In  many  cases,  the  existing  meshes 
no  longer  have  their  original  CAD  geometry  file,  so  it  is  impossible  to  generate  high- 
order  meshes  even  if  software  to  do  so  is  available.  To  address  these  issues,  the  general- 
purpose  low-to-high-order  mesh  converter  meshCurve  was  created.  The  software  is  a 
cross-platform,  GUI-based  tool  to  easily  and  efficiently  perform  low-to-high-order  mesh 
conversion.  Access  to  the  CAD  geometry  is  unnecessary  because  the  software  uses  the 
mesh  itself  to  infer  the  surface  curvature,  while  automatically  identifing  and  preserving 
feature  curves.  Since  high-order  meshing  is  an  issue  for  multiple  simulation  disciplines, 
we  anticipate  that  meshCurve  may  also  be  useful  in  disciplines  beyond  CFD. 

1.2  Program  Objectives 

The  objective  of  this  work  is  to  develop  a  robust,  accurate  and  efficient  high-order  CFD 
tool  for  real  world  engineering  problems.  In  particular,  the  specific  aims  of  this  research 
are  identified  as  follows: 

•  Develop  a  dual-consistent  discrete  adjoint  equation  for  the  CPR  method. 
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Develop  an  efficient  and  accurate  automated  adjoint-based  error  estimation  method 
for  the  CPR  formulation. 

Implement  an  adjoint-based  adaptive  solver  with  the  CPR  method,  which  can  han¬ 
dle  arbitrary  discretization  orders,  and  compare  the  different  adaptation  strategies. 

Demonstrate  the  importance  of  hp-adaptation  for  the  high-order  CPR  method  to 
aerodynamic  flows  and  engineering  problems. 

Evaluate  SGS  models  with  3D  turbulent  flows  using  FR/CPR  method. 

A  priori  and  a  posteriori  evaluation  of  SGS  models  with  Burgers’  Equation  and 
Euler  Equations  using  FR/CPR  method. 

Demonstrate  the  accuracy  and  efficiency  of  LES  with  the  high-order  CPR  method 
to  aerodynamic  flows  and  apply  it  to  a  wide  range  of  engineering  applications. 

Develop  a  GUI-based,  user-friendly,  linear  to  high-order  mesh  converter. 
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2  Results  and  Discussion  on  Adjoint-based  Adaptive 
High-order  Differential  Formulation 

2.1  Review  of  the  Correction  Procedure  via  Reconstruction  Method 

The  CPR  method  [54,  108]  can  be  derived  by  transforming  a  weighted  residual  form  into 
a  differential  one.  Consider  a  hyperbolic  conservation  law 

^+V-F(Q)=  0,  (1) 

with  proper  initial  and  boundary  conditions,  where  Q  is  the  state  vector,  and  F  =  (F,  G) 
is  the  flux  vector.  Assume  that  the  computational  domain  £1  is  discretized  into  N  non¬ 
overlapping  triangular  elements  {V,}fL1.  Let  W  be  an  arbitrary  weighting  or  test  function. 
The  weighted  residual  formulation  of  Eq.  1  on  element  V,  can  be  expressed  as 

dQ  =  0-  (2) 

Let  Qi  be  an  approximate  solution  to  the  analytical  solution  Q  on  V,.  On  each  element, 
the  solution  belongs  to  the  space  of  polynomials  of  degree  k  or  less,  i.e.,  Qi  £  Pk(Vi). 
After  applying  integration  by  parts  twice  to  the  flux  divergence  and  replacing  the  normal 
flux  term  with  a  common  Riemann  flux  F™om  in  the  above  equation,  we  get 

J  ^IVdb!  +  fwV-  F(Qi) dfi  +  J  W  [F™om  -  Fn(Qi)\  dS  =  0.  (3) 

Here  the  common  Riemann  flux  F™orn  is  defined  as 

F?orn  =  F?om(Qi,Qi+,n),  (4) 

where  Qi+  denotes  the  solution  outside  the  current  element  Vi,  and  the  normal  flux 
Fn(Qi)  at  the  interface  is 


Fn{Qi)  =  F(Qi)  ■  ft. 

In  order  to  eliminate  the  test  function,  the  boundary  integral  above  is  cast  as  a  volume 
integral  via  the  introduction  of  a  “correction  field”  on  V,  Si  £  Pk(Vi), 


where  [Fn]  =  F?om 
we  obtain 


I  Vi 


W5i&Q,=  /  W[Fn]dS, 


iavi 


(5) 


Fn(Qi)  is  the  normal  flux  difference.  Substituting  Eq.5  into  Eq.3, 


fv  (j^  +  v  F(Qi)  +  ^w<in  =  o.  (6) 

If  the  flux  vector  is  a  linear  function  of  the  state  variable,  then  V  •  F(Qi)  £  Pk.  In 
this  case,  the  terms  inside  the  square  bracket  are  all  elements  of  Pk.  Because  the  test 
space  is  selected  to  ensure  a  unique  solution,  Eq.6  is  equivalent  to 

+  V  •  ?(Qi)  +*i  =  0-  (7) 

at 

For  nonlinear  conservation  laws,  V  •  F(Qi)  is  usually  not  an  element  of  Pk.  As  a 
result,  Eq.6  cannot  be  reduced  to  Eq.7.  In  this  case,  the  most  obviously  choice  is  to 
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project  V  •  F(Qi)  into  Pk .  Denote  Il(V  •  F(Qi ))  as  a  projection  of  V  •  F(Qi)  to  Pk .  One 
choice  is 


/  n(V  •  F(Qi))Wdfl  =  [  V  •  F(Qi)W&Q..  (8) 

JVi  JVi 

Then  Eq.6  reduces  to 

^  +  n(V  •  F(Qi))  +Si  =  0.  (9) 

Next,  let  the  DOFs  be  the  solutions  at  a  set  of  solution  points  (SPs)  {rij}  (j  varies 
from  1  to  K  =  {k  +  1  )(k  +  2)/2).  Then  Eq.9  is  true  at  the  SPs,  i.e., 

-  +  nj(V  •  F(Qi))  +  Sij  =  0,  (10) 

where  n,  (V  •  F(Q;))  denotes  the  values  of  Il(V  •  F(Qi ))  at  SP  j.  The  efficiency  of  the 
CPR  approach  hinges  on  how  the  correction  field  <5,  and  the  projection  Il(V  •  F(Qi)) 
are  computed.  Two  approaches,  the  Lagrange  polynomial  (LP)  and  the  chain-rule  (CR) 
formulations,  were  suggested  to  compute  the  projection  of  the  flux  divergence  [108].  For 
the  LP  approach,  the  flux  is  assumed  to  belong  to  the  polynomial  space  of  degree  k,  e.g. 
F,  G  £  Pk.  Therefore,  the  projection  of  the  flux  divergence  is  expressed  as 

n(v-F)  =  V  -  {Y^LjFj), 

j 

where  Lj  is  the  Lagrange  interpolation  polynomial  defined  on  SP  j .  For  the  CR  approach, 
the  flux  divergence  V  •  F  is  assumed  to  belong  to  the  polynomial  space  Pk,  which  can  be 
expressed  as 


n(V-F)=^Lj(—  -VQ),-. 

j 

Note  that,  for  a  linear  conservation  law,  the  LP  and  CR  approaches  lead  to  the  same 
formulation.  However,  for  a  nonlinear  flux  equation,  the  CR  approach  can  reduce  the 
aliasing  errors  [108]. 

To  compute  5i,  we  define  k  +  1  points  named  flux  points  (FPs)  along  each  interface, 
where  the  normal  flux  differences  are  computed.  We  approximate  (for  nonlinear  conserva¬ 
tion  laws)  the  normal  flux  difference  [F’"]  with  a  degree  k  interpolation  polynomial  along 
each  interface, 


[n/«i^/  =  E[F"]/>^fp>  (u) 

i 

where  /  is  a  face  (or  edge  in  2D)  index,  and  l  is  the  FP  index,  and  Lfp  is  the  Lagrange 
interpolation  polynomial  based  on  the  FPs  in  a  local  interface  coordinate.  For  linear 
triangles  with  straight  edges,  once  the  solution  points  and  flux  points  are  chosen,  the 
correction  at  the  SPs  can  be  written  as 

5i’j  =  jv\  E  ^2,ajjAFn)f,isP  (l2) 

1  l'  f&dVi  l 

where  OLjj,i  are  lifting  constants  independent  of  the  solution  variables,  Sf  is  the  face  area, 
|  Vi  |  is  the  volume  of  V).  The  details  of  how  to  compute  the  lifting  Constantsa^/,;  can  be 
found  in  Ref.  [108]. 
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In  ID,  a  continuous  flux  polynomial  F  which  is  equal  to  the  common  flux  at  the  inter¬ 
faces  can  be  reconstructed  using  a  piece- wise  analytic  flux  polynomial  Fi(x)  =  II (F(Qi)) 
and  a  correction  term  <Ji(x)  as 


Fi(x)  =  Fi(x)  +  <Tj(  x). 

According  to  Ref. [54],  <Ti(x)  should  approximate  the  zero  function  and  satisfy  the  following 
equation 


at(x)  =  [FtT/2  -  Fi(xi-1/2)]gL(x)  +  [F%™2  -  Fi(xi+1/2)]gR(x). 

Here,  gR{x)  and  gR(x)  are  both  degree  k+1  polynomials  called  correction  functions  with 
the  properties 


9L(xi- 1/2)  =  1,  9L{xi+ 1/2)  =  0,  gR(xi- 1/2)  =0,  gR(xi+ 1/2)  =  1- 
Then  the  CPR  method  for  the  ID  conservation  law  can  be  expressed  as 


d-% f  +  +  [FZ?,2  ~  Fiixi-iMixij)  +  [F%?/2  -  Fix^W^j)  =  0. 

A  series  of  correction  functions  with  different  accuracy  and  stability  properties  were 
developed  in  Ref.  [54].  If  the  correction  function  g  is  chosen  as  right  Radau  polynomials, 
the  DG  method  is  recovered  from  the  CPR  scheme.  In  this  case,  the  correction  function 
denoted  by  g^Q  or  gi  is  perpendicular  to  the  degree  k-1  polynomial  space.  Similarly,  a  g2 
correction  function  is  defined,  which  is  perpendicular  to  the  degree  k-2  polynomial  space. 
In  summary,  for  any  integer  m  Js  1,  a  gm  correction  function  can  be  defined  which  is 
perpendicular  to  pk~m .  For  the  sake  of  simplicity,  the  projection  operator  n  is  omitted 
in  the  rest  of  the  paper. 


2.2  The  Dual-Consistent  CPR  Formulation 

Aircraft  design  engineers  are  usually  interested  in  scalar  engineering  outputs  such  as  lift 
or  drag  coefficients  generated  from  CFD  simulations.  An  adjoint  solution  can  directly 
relate  the  local  residual  to  the  engineering  output.  Adjoint  has  been  used  in  a  wide 
range  of  applications  including  optimal  controls,  design  optimization,  data  assimilation 
and  error  estimation.  There  are  two  approaches  to  obtain  an  adjoint  solution.  One 
can  solve  the  continuous  adjoint  equation  which  is  a  partial  differential  equation  using 
any  numerical  method,  or  directly  solve  the  discrete  adjoint  equation  derived  from  the 
discretized  primal  equation.  As  for  the  primal  problem,  a  numerical  scheme  is  defined  as 
a  consistent  method  if  its  discrete  operator  converges  to  the  continuous  operator,  or  the 
exact  solution  satisfies  the  discrete  numerical  formulation  as  the  mesh  size  approaches 
zero.  Similarly,  for  a  dual-consistent  adjoint  formulation,  the  exact  adjoint  solution  from 
the  continuous  adjoint  equation  should  satisfy  the  discrete  adjoint  equation  in  the  limit  of 
vanishing  mesh  size.  The  dual-consistency  of  a  discrete  adjoint  operator  from  a  numerical 
discretization  is  a  key  component  to  ensure  that  the  optimal  convergent  rate  is  achieved 
for  an  engineering  output.  To  establish  a  robust  and  accurate  functional  error  estimation 
procedure  for  the  CPR  method,  a  dual-consistent  CPR  formulation  is  developed  in  this 
section. 


2.2.1  The  Dual  Problem  and  Adjoint-based  Error  Estimation 

Linear  PDEs  First,  we  review  the  dual  problem  in  the  theory  of  output-based  error 
estimation.  A  detailed  discussion  can  be  found  in  a  series  of  articles  [41,  42]  and  references 
therein.  Consider  a  linear  differential  equation 


DISTRIBUTION  A:  Distributipi])  approved  for  public  release 


LQ  =  f  on  fl, 

with  a  homogeneous  boundary  condition,  where  Q  £  V  is  the  solution,  V  is  the  infinite 
dimensional  solution  space  and  /  £  L2(fl).  Suppose  the  output  functional  of  interest  J 
is  given  as  an  inner  product  of  a  smooth  function  g  and  the  solution  Q 

J(Q)  =  (g,Q)  =  [  gQdV, , 

J  o 

over  the  entire  domain  Q.  The  dual  problem  is  introduced  by  adding  a  weighted  residual 
to  the  functional 


J ( Q )  =  (5,  Q)  +  Wh  /  -  LQ). 

where  ip  is  an  arbitrary  function  for  now.  If  the  solution  Q  satisfies  the  linear  differential 
equation,  this  weighted  residual  does  not  affect  the  value  of  the  original  functional.  Denote 
L*  the  adjoint  differential  operator  with  respect  to  L  defined  according  to  ( ip,LQ )  = 
(L*if>,Q)  with  homogeneous  boundary  conditions.  Then  we  have 


J(Q)  =  (g,Q)+WJ)-(<P,LQ)  =  {g,Q)+(ip,f)-(L*ip,Q)  =  (V>, /)-(£>-<?, Q).  (13) 

Let  1 p  be  the  adjoint  solution  computed  using 

L*il>  =  g  onfl. 

Then  the  last  term  in  Eq.  13  vanishes,  and  we  obtain 

J(Q)  =  Wh/)- 

Obviously  J  is  independent  of  Q,  but  depends  on  if)  now.  So  we  call  J  a  function  of  ip, 
i.e., 


J(iP)  =  (iP,f)-(L*iP-g,Q).  (14) 

Therefore,  the  duality  of  the  functional  is 

J(Q)  =  (g,Q)  Q  LQ  =  f  on  fl 

or 

J(ip)  =  (ip,f)  if  L*  ip  =  g  on  ft. 

Suppose  Qh  and  iph  are  the  discrete  primal  solution  and  the  discrete  adjoint  solution 
obtained  with  a  numerical  method,  and  both  of  them  belong  to  the  discrete  solution  space 
Vh ■  Then,  the  discrete  source  term  for  the  primal  equation  is  fh  =  LQ The  functional 
error  can  be  estimated  by 


=  (g,  Qh)  -  (5,  Q)  =  ( g ,  Qh  -  Q) 

=  {L*ip,  Qh  —  Q) 

=  (ip ,L{Qh  -  Q)) 

=  ( iphjh  -  /)  +(ip  -  iphjh  ~  /),  (15) 

computable  error  remaining  error 
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where  the  linearity  of  the  inner  product  and  the  adjoint  definition  are  used.  Note  that 
the  first  term  on  the  RHS  of  Eq.  15  is  defined  as  the  computable  error,  since  it  does 
not  involve  any  analytical  primal  solution.  If  a  numerical  method  possesses  the  Galerkin 
orthogonality  property, 


(vh,  fh~  f)  =  0,  for  all  vh  G  Vh , 

the  computable  error  vanishes  since  we  take  ifh  in  space  V), .  Therefore,  there  is  no  need 
to  evaluate  the  computable  error,  and  the  order  of  the  output  error  only  depends  on  the 
remaining  error  term.  If  we  use  a  degree  k  polynomial  to  approximate  the  primal  and 
adjoint  solution  and  assume  the  optimal  order  of  accuracy  can  be  achieved,  we  obtain 

WQh  -  Oil  =  0{hk+1),  \\ifh  -  if\\  =  0{hk+1). 

Here 1 1 . 1 1  is  a  L2  norm.  In  addition,  for  a  nth  order  differential  PDE,  we  have  \\fh  —  /||  = 
0{hk+1~n).  So  the  order  of  the  computable  error  is  0(h2k+2~n),  which  leads  to  a  su- 
perconvergent  functional  of  order  0{h2k+2~n).  To  estimate  the  output  error,  we  need  to 
evaluate  the  computable  error.  A  common  approach  is  to  use  an  adjoint  solution  from  a 
finer  approximation  space,  e.g.  p  =  k  +  1.  Then  the  computable  error  is  not  equal  to  0, 
and  we  can  use  this  error  estimate  to  correct  the  original  output 

fJcorr  —  ffiQh)  +  {if h j  fh  ~  /)• 

The  convergence  rate  of  this  corrected  output  is  two  orders  higher  than  the  original 
output,  which  is  0{h2k+4~n).  A  rigorous  proof  can  be  found  in  Ref.  [52], 

Nonlinear  PDEs  Consider  a  non-linear  differential  equation 

J\f(Q)  =  0,  onfl.  (16) 

Suppose  an  output  functional  of  interest  is  given  as  ff{Q).  A  dual  problem  is  introduced 
by  defining  a  Lagrangian  of  the  output  with  the  constraint  of  the  solution  Q  satisfying 
the  primal  equation  AT (Q)  =  0 


C  =  J(Q)  +  [  ifAf(Q)d£l.  (17) 

Jo. 

Here  if  €  V  has  two  roles.  First,  if  is  the  adjoint  solution.  Second,  it  also  serves  as  a 
Lagrangian  multiplier.  After  performing  the  linearization  and  enforcing  stationary  of  C 
to  a  solution  perturbation  SQ  G  V,  we  obtain 

SC  =  C'[Q](6Q)  =  J'[Q](5Q)  +  f  ifAf'[Q](6Q)dfl  =  0  V5Q  G  V,  (18) 

■In 

where  the  primed  notation  denotes  Frechet  linearization  with  respect  to  an  argument  in 
the  square  bracket.  Eq.  18  defines  the  dual  problem  in  a  variational  form  by  finding  if 
such  that 


J'[Q\(v)+  [  ifAf'[Q](v) dH  =  0  VuGF.  (19) 

Jn 

Let  Qh  denotes  an  approximate  solution  to  the  analytical  solution  Q.  The  difference 
between  them  can  be  interpreted  as  a  solution  perturbation  Qh  =  Q  +  SQ.  The  output 
error  defined  as  SJ  =  J{Qh)  —  J{Q)  can  be  estimated  by  the  adjoint  weighted  residual 
method 
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SJ  «  J'[Q\(SQ) 


*l>Af'[Q]{6Q)<m  »  -  /  V’(A^(Q?l)-AA(Q))dri  =  - 


ipAf{Qh)dfl. 

(20) 


2.2.2  The  Continuous  Adjoint  Equation 

We  consider  the  following  conservation  law 

A/"(Q)  =  V  •  ?{Q),  (21) 

as  an  example  to  develop  the  continuous  adjoint  equation.  Eq.  19  leads  to 

J,[Q](u)  +  ^V’(^  (V-/))  (u)dO  =  0  VreV.  (22) 

Suppose  the  output  functional  J(Q)  consists  of  surface  (dQ)  and  volume  (0)  integrals 


J{Q)=  [  gQdtt+  [  jT(Q)ds , 

J  n  Jan 

where  jT  is  a  boundary  operator.  Substituting  the  definition  of  J  into  Eq.  22,  and 
performing  integration  by  parts,  we  get 


[{9~^'%)dn+  f  (§§+^M'^)ds  =  0' 

.In  dQ  Jan  dQ  oQ 


VV’- 


df 

dQ 


=  9, 


(23) 


which  is  a  linear  partial  differential  equation  for  ip,  and  the  corresponding  boundary 
conditions  are  from  the  surface  integral 


djT 

dQ 


■  n)ds  =  0. 


(24) 


2.2.3  The  Dual-Consistent  CPR  Formulation 

The  discrete  adjoint  equation  is  obtained  directly  by  linearizing  the  discretized  primal 
equation.  Consider  a  discretized  formulation  of  the  primal  Eq.  16 

A fh{Qh)  =  0,  Qh  €  Vh 

with  mesh  h,  and  a  discrete  solution  perturbation  SQ h  £  V), .  Linearizing  the  discrete 
residual  7 Zh  and  the  discrete  output  Jh,  we  get 

SKh  =  nh(Qh  +  SQh )  —  lZh{Qh)  ~  -wTT-^Qh 

SJh  =  Jh(Qh  +  SQh)  -  JhiQh)  ~  -J^r^Qh- 

oQh 

The  discrete  adjoint  % ph  is  defined  as  the  sensitivity  of  output  perturbation  SJh  to  the 
primal  residual  perturbation  SIZh 

SJh  =  -ffiSKh, 
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forms 


Figure  1:  The  Cl  adjoint  for  a  subsonic  NACA0012  airfoil 


which  leads  to 


f ^SQh  «  SJh  =  «  -t’h^SQh- 

oQu  oQh 

We  obtain  the  discrete  adjoint  equation  by  canceling  the  solution  perturbation  SQh 


dnh 


4>h  = 


djh 


(25) 


dQh  dQh 

For  numerical  methods  with  a  weak  form  such  as  FEM  or  DG  methods,  after  choosing 
a  proper  basis,  Eq.  25  is  equivalent  to  its  variational  formulation.  A  detailed  derivations 
can  be  found  in  Ref.  [47,  30].  So  the  discrete  adjoint  equation  for  numerical  methods  in 
semi-linear  form  is  consistent  with  the  continuous  adjoint  equation.  However,  this  is  not 
true  for  numerical  methods  in  a  differential  form  such  as  the  CPR  method,  which  does 
not  possess  a  variational  form. 

Substituting  a  pointwise  residual  ry.j  defined  on  each  solution  point  j  of  cell  i  arising 
from  a  differential  scheme,  we  obtain 


VV^-  -  =  — 

4-  4^  dQl  3  dQt  ’ 

*  3 


(26) 


where  l  is  the  global  index  of  the  DOFs  in  the  entire  domain.  Figure  la  shows  the  x- 
momentum  component  of  the  lift  adjoint  for  a  subsonic  airfoil  using  the  discrete  adjoint 
equation  with  the  CPR  method.  It  has  a  very  oscillatory  distribution  in  each  cell,  which 
indicates  dual-consistency  violations. 

Since  the  CPR  method  is  not  in  a  variational  form,  its  discrete  adjoint  equation 
should  be  directly  derived  from  the  linearized  Lagrangian.  Assume  the  adjoint  solution 
belongs  to  the  same  space  of  the  primal  solution,  the  adjoint  variable  t/>j  of  cell  i  can  be 
approximated  using  the  Lagrange  basis  Lj 


i’i  =  5Z  Li^r 

3 

Directly  discretizing  the  linearized  Lagrangian,  Eq.  19,  with  a  quadrature  rule,  we  obtain 


-EE 

*  3 


dr. 


dQi 


dj 
dQi ’ 


(27) 
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where  u3  and  \Jtj\  are  the  quadrature  weight  and  the  element  Jacobian  at  the  solution 
point.  Comparing  with  Eq.  26,  the  following  relation  can  be  derived  from  the  discrete 
adjoint  ipij  and  the  discrete  adjoint  ipij  in  the  integral  form 

It  is  obvious  that  the  discrete  adjoint  formula  for  a  numerical  scheme  in  a  differential  form 
is  not  consistent  with  the  continuous  adjoint  equation.  The  inconsistent  adjoint  is  related 
to  the  consistent  counterpart  in  terms  of  the  quadrature  weight  u>  and  the  cell  Jacobian 
|  J\.  This  integral  equation  can  be  interpreted  as  an  explicitly  defined  variational  form  for 
the  CPR  method.  In  this  paper,  we  call  it  the  dual  consistent  discrete  adjoint  formula. 
Figure  lb  shows  the  consistent  discrete  adjoint  solution  with  the  CPR  method.  Clearly 
the  consistent  adjoint  solution  is  smooth. 

2.2.4  Analysis  of  Dual-Consistency  for  the  CPR  Method 

For  a  dual  consistent  discretization,  the  discrete  adjoint  equation  should  approach  the 
continuous  adjoint  partial  differential  equation  when  the  mesh  size  diminishes.  In  other 
words,  the  analytic  primal  solution  Q  and  the  analytic  dual  solution  ip  should  satisfy 
the  discrete  adjoint  equation  when  the  mesh  size  goes  to  zero.  Substituting  the  discrete 
residual  of  the  CPR  method  into  the  discrete  adjoint  equation  for  the  conservation  law 
leads  to 


3  /  \  3  7 

E£^(v  ■  F(Q)ij  +  Sijj  Uj\Ji,j\'4>i,j  —  3Ql  ’  (^8) 

*  3 

where  Sl  :l  is  the  correction  term.  For  a  ID  CPR  formulation,  the  above  equation  can  be 
further  expressed  as 


dFi,k  dLk(£j) 

~dQi  df~ 


dQi 


[F]i+igR(£,j)  J  i,j 


dj 

dQi 


(29) 

Here,  g'L,  g'R  and  t/>,  belong  to  Pfc(fl,;).  ^  is  a  degree  k  —  1  polynomial.  Therefore,  the 
degree  of  the  integrand  is  at  least  2k.  Assume  that  the  quadrature  rule  defined  on  the 
solution  points  is  exact  at  least  for  a  degree  2k  polynomial.  Then  we  have 


?  /„.  *  ?  W^X+^L  *  (s||F1-^  +  jIjI'W*)  «•*>  =  - 


3J 

dQi 


Performing  integration  by  parts  on  the  LHS,  we  obtain 


LHS  =  -J_ l 


f  dF  dip  ^  dF  ]i+! 

/0,ae<fcda:  +  5>i 


^  ,  .  a 


V’igLdC 


A-  d^+h  I  ,  ,i 


V’igRdC 


(30) 


Recall  that,  for  the  CPR  method,  the  correction  functions  satisfy 
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9l(~  !)  =  9r(  1)  =  1 
9l{  1)  =  9r{~  !)  =  0. 


Furthermore,  a  DG  correction  function  gdg  of  degree  k  +  1  is  perpendicular  to  Pk  1 ,  i.e. , 

f  9dg(Om<%  =  0,  .VtfeP*"1. 

Therefore,  if  a  DG  correction  function  is  used,  Eq.  30  can  be  simplified  to 


LHS  =  — 


Z^ 


dF  dip  i  dF  ,i+  i 

dQtedX  +  ^dQ  l*-l 


-E 


9Q 


^-i+E 


5Q 


V’i+I,.,  (31) 


or  equivalently, 


lhs  =  -J2 


dF  dip 
dQ  dx 


dx  ■ 


Z^ 


f)  rpcom 
dti+ 1 

dQ 


dFc 


V’i+i  - 


z^ 


dQ 


-ipi-l 


(32) 


The  first  term  is  the  governing  equation  for  the  continuous  adjoint  equation,  which  van¬ 
ishes  for  the  analytic  adjoint  solution.  Notice  that 


QP'com  Qp'com 


Z^ 


f)  TFcom 
N+i 

dQ 

f)  r?com 
N+  i 


QP'com  N—  1  Qjpcom 

V’at+i — +  Ew^1)  -ifc+i (-1))-  i+5 

i=l 


V’AT+i  — 


9Q 

aPT" 


5Q 


-ip  i 


dQ  rjv+5  3Q 

when  ipi(l)  =  ipi+i ( — 1)  for  the  analytical  adjoint  solution.  So  the  final  equation  is 


jj'com  f)  jj'com 


dj 
dQ ' 


(33) 


As  discussed  in  Ref.  [47,  75],  a  well-defined  dual-consistent  boundary  flux  should  only  be 
a  function  of  the  boundary  state  Fcom  =  Fcom(Qbc(QiJ))1  and  a  properly-defined  dual- 
consistent  output  functional  leads  to  Eq.  33  and  is  consistent  with  the  dual  boundary 
condition  for  the  continuous  adjoint  equation  (Eq.  24) .  Therefore,  Eq.  28  can  be  satisfied 
exactly  with  an  analytical  adjoint  solution.  A  similar  procedure  can  be  applied  to  the 
system  of  equations  in  2D.  The  corresponding  equation  to  Eq.  33  is 


ip 


gpcom 

ixrds= 


dj 
dQ ’ 


(34) 


which  will  be  used  to  analyze  the  dual  consistency  for  the  2D  linear  wave  equation  and 
the  Euler  equations  in  the  next  section. 

Based  on  the  analysis,  the  key  factors  to  ensure  a  dual-consistent  CPR  formulation 
are  summarized  next. 


1.  In  order  to  ensure  the  integral  accuracy  of  the  discrete  adjoint  equation  in  a  vari¬ 
ational  form,  an  accurate  quadrature  rule  defined  on  the  solution  points  should  be 
exact  for  a  degree  2k  polynomial.  Recall  that  a  k+1  point  Gaussian  quadrature  rule 
can  yield  an  exact  integration  of  a  degree  2k  +  1  polynomial.  Therefore,  the  Gauss 
quadrature  points  are  preferred  as  the  solution  points.  The  Lobatto  quadrature  rule 
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can  only  integrate  a  degree  2k  —  1  polynomial  exactly.  If  the  Lobatto  points  are 
used  as  the  solution  points,  the  CPR  formulation  will  have  an  accuracy  loss  for  the 
discrete  adjoint  solution  and  the  corresponding  error  estimation. 

2.  The  correction  function  g  in  the  CPR  method  must  be  perpendicular  to  the  deriva¬ 
tives  of  the  adjoint  solution 


I  gip'dCl  =  0. 
■In 


(35) 


If  we  assume  that  the  discrete  adjoint  ip  belongs  to  the  same  space  of  the  primal 
solution,  the  degree  of  the  derivative  ip'  is  k  —  1.  Then  the  only  qualified  correction 
function  is  the  DG  correction  function  gdg,  which  is  perpendicular  to  Pk _1.  On 
the  other  hand,  the  degree  of  the  adjoint  solution  ip  is  determined  by  the  specific 
correction  function.  Suppose  we  use  the  g2  correction  function,  which  is  only  per¬ 
pendicular  to  Pk~2.  To  satisfy  this  condition,  the  degree  of  the  discrete  adjoint 
solution  ip  is  automatically  degenerated  to  k  —  1.  In  summary,  for  a  CPR  scheme 
with  a  correction  function  gm,  the  discrete  adjoint  ip  G  pk+l~m. 

3.  The  Lagrange  polynomial  (LP)  approach  is  required  to  evaluate  the  flux  divergence 
term  in  the  CPR  formulation,  instead  of  the  chain  rule  (CR)  approach.  This  re¬ 
quirement  is  very  similar  to  the  conservation  requirements  for  the  CPR  scheme. 
Therefore,  a  similar  fix  can  be  obtained  by  following  the  conservation  fix  of  the  CR 
approach  [35].  However,  the  numerical  results  in  the  next  section  show  that  the 
dual-inconsistent  violation  by  the  CR  approach  is  relatively  weak. 

4.  A  properly  defined  common  numerical  flux  on  the  boundaries,  and  a  well-defined 
output  functional,  are  critical  in  making  the  boundary  terms  in  the  adjoint  equation 
vanish. 


2.2.5  Numerical  Tests 

Linear  Advection  Equation  First,  the  dual-consistent  CPR  discretization  and  the 
error  estimation  are  verified  with  a  first-order  hyperbolic  partial  differential  equation 
(PDE).  Consider  a  2D  linear  wave  equation 


V  •  cQ  =  /,  itil 
with  a  Dirichlet  boundary  condition 

Q(x)  =  B{x),  x  €  dfl~ 

given  on  the  inflow  boundaries  dQ~  ={it  <9fl|  n-  c  <  0},  where  n  is  the  outward  surface 
normal  direction  and  c  is  the  advection  velocity  (Figure  2a)  prescribed  as 

c  —  (e~x,  ey). 

The  solution  is  set  to  be 


7t(gX  —  l) 

Q(ff=e"am(  1  .  j), 

e  —  1 

which  determines  the  source  term  f  by  the  method  of  manufactured  solutions  [50] . 

The  simulation  is  performed  on  a  unit  square  ff  =  [0,1] 2  filled  with  quadrilateral 
elements.  The  output  of  interest  is  defined  on  the  outflow  boundaries  <9fl+  =  {x  £ 
<9fl|  n  ■  c  >  0}  as 
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where  the  weighting  function 


J{Q)=  [  g(x)(n  ■  cQ)ds, 
JdCl+ 


a{x) 


The  exact  value  of  the  output  is  J  =  2e. 
In  this  case,  Eq.  34  leads  to 


7T  ex~y 
e  —  1 


ion+ 


Qjpcom  r 

ip  a<n  ds  =  /  g(x)(n- cQ)ds, 

dQ  oQ  Jqq+ 


which  is  equivalent  to  the  boundary  condition  for  the  continuous  adjoint  equation 


/  (n  ■  c)  ip  +  (n  ■  c)g(x)ds  =  0. 

JdU+ 

The  common  flux  Fcom(QBc(QL),n)  =  n • cQl  is  chosen  for  the  boundaries  to  ensure 
that  the  CPR  method  for  the  linear  wave  equation  with  this  output  is  a  dual-consistent 
formulation. 

The  contours  of  the  primal  solution  and  the  adjoint  solution  are  shown  in  Figure  2. 
Based  on  the  dual-consistency  analysis,  we  choose  the  Gauss  points  as  the  solution  points 
and  the  DG  correction  function  gdg .  Figure  3  shows  the  convergent  rate  of  the  primal 
solution,  the  output  functional  and  the  error  estimate  using  the  CPR  method.  The  opti¬ 
mal  order  0(hk+1)  is  obtained  for  the  primal  solution.  Super  convergence  of  the  output 
functional  and  the  error  estimate  is  good  indicator  of  a  dual-consistent  discretization. 
For  this  case,  a  super  convergence  of  0(h2k+1)  is  observed  for  both  of  the  output  error 
and  the  adjoint-based  output  error  estimate.  The  corrected  output  converges  2  orders 
faster  than  the  original  output,  which  is  0(h2k+3).  The  results  of  the  CPR  method  with 
Lobatto  points  and  gdg  is  shown  in  Figure  4a.  Comparing  with  the  results  of  the  Gauss 
points,  the  accuracy  loss  of  the  quadrature  rules  defined  on  the  Lobatto  points  leads  to 
one  order  loss  in  the  output  functional,  the  error  estimate  and  the  corrected  functional. 

Now  we  test  the  influence  of  the  correction  functions  on  the  output  functional  and 
the  error  estimate.  As  discussed  in  [54],  for  a  integral  m  >  1,  a  correction  function  gm 
of  degree  k  +  1  is  perpendicular  to  Pfc_m,  and  the  Fourier  analysis  indicates  that  the 
order  of  the  corresponding  scheme  is  0{h2k+2~m).  The  previous  analysis  shows  that  the 
approximation  order  of  the  adjoint  solution  ip  is  determined  by  the  correction  function 
gm  as  0(hk+1~m). 

Figure  4  shows  the  convergence  rates  using  the  gi  and  g%  correction  functions  with  the 
Gauss  points  as  the  solution  points.  Similar  relationships  are  obtained  between  the  order 
of  accuracy  of  the  functional  outputs,  the  error  estimate  and  the  correction  functions. 
For  the  advection  equation,  the  output  functional  and  the  adjoint-based  error  estimate 
with  the  Gauss  points  and  gm  are  accurate  to  order  0(h2k+2~m).  Furthermore,  the 
corrected  outputs  of  the  corresponding  schemes  are  accurate  to  order  0(h2k+4~m).  The 
convergence  rates  with  the  different  CPR  schemes  are  summarized  in  Table  1.  Those 
results  are  consistent  with  the  analysis  in  section  2.2.3. 

For  a  linear  partial  differential  equation,  the  discrete  adjoint  based  output  error  esti¬ 
mation  should  recover  the  exact  output  error  obtained  by  a  finite  difference  method.  Table 
2  compare  the  true  output  error  between  k  =  1  and  =  =  2  results  and  the  adjoint-based 
output  error.  The  difference  between  them  are  within  machine  zero. 
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(a)  Primal  solution 


(b)  Adjoint  solution 


Figure  2:  The  primal  and  adjoint  solution  of  the  linear  advection  equation 


SPs 

9 

primal  L2  err. 

output  error 

output  error  est. 

corrected  output  error 

Gauss 

DG 

k+1 

2k+l 

2k+l 

2k+3 

Lobatto 

DG 

k+1 

2k 

2k 

2k+2 

Gauss 

G2 

k+1 

2k 

2k 

2k+2 

Gauss 

G3 

k+1 

2k- 1 

2k- 1 

2k+l 

Table  1:  Convergence  rates  for  the  linear  wave  equation  with  different  schemes 


The  Supersonic  Vortex  Transportation  Problem  Second,  we  test  the  present 
adjoint-based  error  estimation  for  the  Euler  equations  with  curved  elements.  The  problem 
we  consider  is  a  2D  supersonic  vortex  transported  in  a  circular  sector.  The  computational 
domain  is  defined  on  a  section  of  an  annulus  with  the  inner  radius  of  Tin  =  1  and  the 
outer  radius  of  1.384.  The  initial  mesh  consists  of  k  =  4  quadrilateral  elements  and  is 
shown  in  Figure  5a.  The  isentropic  vortex  rotates  around  the  center  of  the  circular  sector. 
The  density  p  is  only  a  function  of  the  radius  r  (Figure  5b) 

P(r)  =  PU 1  +  -  r-f))W, 

2  rz 

where  7  is  the  ratio  of  heat  capacities  and  the  remaining  parameters  are  the  flow  conditions 
on  the  inner  surface  chosen  as  pin  =  2,  Min  =  2  and  pin  =  .  The  other  variables  can 

be  computed  with  the  isentropic  relations.  Characteristic  boundary  conditions  with  the 
analytical  solution  are  used  at  both  the  inlet  and  the  outlet,  and  slip  wall  boundary 
conditions  are  applied  on  the  inner  and  the  outer  boundaries.  The  output  of  interest  is 
the  force  in  the  x  direction  on  the  inner  surface,  where  the  pressure  is  equal  to  so  the 
exact  value  of  the  output  is  J  =  —  ^ . 

The  output  boundary  operator  jT  is  defined  as  jT  =  pn  ■  idir ,  where  idir  =  [1,0]. 
Therefore,  the  corresponding  equation  to  Eq.  33  in  2D  is 

f  TdFc°m  d  f  dp  v  w 

/  V  O/O  ds  =  ~-^n  /  -K^{n-tdir)ds, 

id n  9Q  9Q  Jon  9Q 

which  is  equivalent  to  the  boundary  condition  for  the  continuous  adjoint  equation 
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(c)  P3  (d)  P4 


Figure  3:  Results  for  the  linear  advection  equation  using  a  CPR  scheme  with  the  Gauss 
points  and  gdg 
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chosen  as  the  common  flux  on  the  boundaries.  This  ensures  that  the  CPR  method  for 
the  Euler  equations  with  this  output  definition  is  a  dual-consistent  formulation. 

The  adjoint  solution  with  the  dual-consistent  boundary  conditions  and  dual-inconsistent 
boundary  conditions  are  shown  in  Figure  5.  The  dual-inconsistent  boundary  conditions 
generated  some  spurious  oscillations  near  the  wall,  while  the  adjoint  solution  from  the 
dual-consistent  boundary  conditions  was  very  smooth.  Figure  6  displays  the  solution  er¬ 
ror,  the  output  error  and  the  adjoint-based  error  estimate.  The  optimal  order  of  accuracy 
k  +  1  in  L2  norm  was  obtained  by  both  the  dual-consistent  and  inconsistent  boundary 
conditions  in  the  primal  solution.  However,  a  super  convergence  of  order  2k  +  1  was 
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Error  Error  Error 


(a)  Lobatto  points  with  gdg  as  the  correction  function 


SQRT(DOFs)  SQRT(DOFs) 

(c)  Gauss  points  with  <73  as  the  correction  function 

Figure  4:  Results  for  the  linear  advection  equation  with  different  correction  functions 
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(c)  Adjoint  solution  V']  with  dual-(d)  Adjoint  solution  ipi  with  dual- 
consistent  BCs  inconsistent  BCs 

Figure  5:  The  primal  and  adjoint  solution  for  the  supersonic  vortex  transportation  prob¬ 
lem  (Gauss  points,  gdg,k  =  3) 


(a)  L 2  density  error 


(b)  Output  error 


Figure  6:  The  results  of  the  supersonic  vortex  transportation  problem  using  the  dual- 
consistent  BC  and  the  dual- inconsistent  BC  (Gauss  points,  gdg) 
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Cells 

Jh(Qh)  ~  Jh{Qh ) 

-WhVRhiQZ) 

4 

-7.85516345138E-002 

-7. 855163451 17E-002 

16 

-1.25394717994E-002 

-1.25394717999E-002 

64 

-1.83796281052E-003 

-1.83796281080E-003 

256 

-2.63400045419E-004 

-2.63400045485E-004 

1024 

-3.58710248100E-005 

-3.58710248149E-005 

Table  2:  The  output  error  by  the  finite  difference  method  and  adjoint-based  error  esti¬ 
mation  (coarse  space  k  =  1,  fine  space  k  =  2) 
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V  •  F 

primal  L2  err. 

output  error 

output  error  est. 

corrected  output  error 

Gauss 

DC 

LP 

k+1 

2k+l 

2k+l 

2k+3 

Gauss 

DG 

CR 

k+1 

2k+l 

2k+l 

2k+3 

Lobatto 

DG 

LP 

k+1 

2k 

2k 

«  2k  +  1 

Lobatto 

DG 

CR 

k+1 

2k 

2k 

2k+l 

Gauss 

G2 

LP 

k+1 

2k 

2k 

2k+2 

Gauss 

G3 

LP 

k+1 

2k- 1 

2k- 1 

2k+l 

Table  3:  Order  of  accuracy  for  the  supersonic  vortex  transportation  problem 


observed  for  the  output  with  the  dual-consistent  boundary  condition  only.  The  spuri¬ 
ous  adjoint  oscillation  caused  by  the  dual-inconsistent  boundary  condition  degraded  the 
adjoint  solution,  and  destroyed  the  super  convergence  property  of  the  output  functional. 

Figure  7  and  8  show  the  convergence  rates  of  the  CPR  method  with  different  solution 
points  using  both  the  LP  and  the  CR  approaches.  Similar  results  were  obtained  with 
the  LP  and  CR  approaches  for  the  Gauss  points:  a  super-convergence  of  order  2k  +  1  for 
the  output  functional  and  the  error  estimate  and  a  super  convergence  of  2k  +  3  for  the 
corrected  output.  This  indicates  that  the  dual-consistency  violation  of  the  CR  approach 
is  relatively  weak,  and  does  not  affect  the  adjoint-based  error  estimate.  However,  with 
the  Lobatto  points  and  the  LP  approach,  accuracy  loss  did  occur.  The  super  convergence 
of  the  corrected  output  is  lost  for  the  k  =  1  scheme.  The  CPR  schemes  with  Lobatto 
points  and  the  CR  approach  can  reduce  the  alias  error  generated  by  the  non-linear  fluxes. 
Even  though  the  CR  approach  is  not  fully  dual-consistent,  the  super  convergence  rates 
are  recovered  for  the  Lobatto  points,  whose  output  functional  and  the  error  estimate  are 
accurate  to  a  super  convergence  order  of  2k  and  the  order  of  the  corrected  output  is 
2 k  +  1.  Figure  9  shows  the  results  of  the  CPR  method  with  different  correction  functions. 
Similar  results  as  the  linear  wave  equation  are  obtained.  For  the  Euler  equations,  the 
output  functional  and  the  adjoint-based  error  estimate  with  the  Gauss  point  and  gm  are 
accurate  to  order  0(h2k+2~m).  Furthermore,  the  corrected  outputs  of  the  corresponding 
schemes  are  accurate  to  order  0(h2k+4~m).  The  convergence  rates  with  the  different 
CPR  schemes  are  summarized  in  Table  3.  The  results  of  this  test  case  indicate  that  the 
dual  consistent  formulation  performs  as  expected  for  a  non-linear  equations  and  curved 
elements. 

Inviscid  Flow  over  the  NACA0012  Airfoil  The  last  test  case  in  this  section  is  a 
subsonic  flow  over  a  NACA0012  airfoil  with  a  free-stream  Mach  number  of  =  0.5 
and  an  angle  of  attack,  a  =  2°.  It  is  used  to  further  assess  the  accuracy  of  the  adjoint- 
based  error  estimation  for  a  problem  with  a  geometric  singularity.  The  output  of  interest 
is  chosen  as  the  lift  or  drag  of  the  airfoil.  The  contours  of  the  Cl  adjoint  and  the  Cd 
adjoint  are  shown  in  Figure  lb  and  Figure  10.  In  this  test  case,  the  error  in  the  functional 
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(b)  Gauss  points,  CR  and  gfig 

Figure  7:  The  results  of  the  supersonic  vortex  transportation  problem  with  the  Gauss 
points 
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(b)  Lobatto  points,  CR  and  g^g 


Figure  8:  The  results  of  the  supersonic  vortex  transportation  problem  with  the  Lobatto 
points 
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Figure  9:  The  results  of  the  supersonic  vortex  transportation  problem 
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Figure  10:  The  3rd  component  of  the  Cd  adjoint  of  an  Inviscid  NACA  0012  Airfoil  at 
Moo  =  0.5,  a  =  2° 


JH(QH)-MQh)  is  computed  using  p-enrichment  from  k  =  2  to  k  =  3  and  the  effectivity 
of  the  error  estimation  is  defined  as 


^ H  Jh(Qh)  ~  Jh{Qh) 

Due  to  the  geometry  singularity,  the  super  convergence  of  the  output  functional  and 
its  error  estimate  is  lost.  Table  4  shows  the  results  with  4  levels  of  uniformly  refined 
meshes.  Note  that  the  error  of  the  initial  estimates  on  the  very  coarse  meshes  is  large; 
however,  the  effectivity  index  rjfj  approaches  unity  as  the  mesh  is  refined  for  both  Cl  and 
Cd-  In  addition,  the  error  estimate  with  the  Cd  adjoint  is  more  accurate  than  with  the 
Cl  adjoint.  This  is  due  to  the  fact  that  the  regularity  of  the  Cl  adjoint  is  low  because  of 
the  singularity  along  the  stagnation  streaming  line  from  the  leading  edge,  while  the  Cd 
adjoint  is  relatively  smooth. 


Cells 
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Jh(Qh)  -  Jh{Qh ) 
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1.011 

4480 
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-1.68E-03 
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2.29E-04 

2.29E-04 

1.002 

17920 

-3.08E-04 

-3.24E-04 

1.049 

2.13E-05 

2.13E-05 

0.999 

Table  4:  Adjoint-based  error  estimate  for  an  inviscid  NACA  0012  airfoil  at  =  0.5,  a  = 
2°  (coarse  space  k  =  2,  fine  space  k  =  3) 


2.3  Adjoint-based  Error  Estimation  and  H-adaptation 

2.3.1  Adjoint-based  H-adaptation 

An  adjoint-based  error  estimation  relates  a  specific  functional  output  error  directly  to  the 
local  residual  error.  Therefore,  it  can  be  used  to  construct  a  very  effective  error  indicator 
to  drive  an  adaptive  procedure  for  any  engineering  output.  From  Eq.  20,  the  output  error 
can  be  estimated  by  performing  a  quadrature  rule  as 

SMQh) «  -  EE  I  (Q h  )’ 

*  3 
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Figure  11:  The  procedure  of  the  adjoint-based  h-adaptation 


The  continuous  adjoint  solution  ip  is  approximated  by  solving  iph  on  the  finner  space 
through  enriching  the  degree  of  the  solution  polynomial.  The  finer  solution  Q /,  is  obtained 
by  performing  several  iterations  of  GMRES  relaxation  after  prolongating  from  the  coarse 
solution  Qh 


Qh  =  Ih  Qh, 

with  an  injection  operator  1^  .  The  adjoint-based  local  error  indicator  r/,:  used  in  the 
present  study  is  defined  by  taking  an  absolution  value  of  the  elemental  output  error 
contribution 


Vi  =  \  (V> h  -  Ihi>H)itj  n,j(Qh)  |. 

3 

Here,  to  achieve  a  better  estimate,  the  adjoint  defect  between  the  coarse  level  and  fine 
level  iph  —  I^ipH  is  used.  For  a  system  of  equations,  the  local  error  indicator  is  formed 
by  summing  together  every  component’s  contribution  to  the  functional  error  estimate. 

Figure  11  shows  the  procedure  of  the  adjoint-based  h-adaptation.  A  fixed-fraction 
hanging-node  h-adaptation  strategy  is  used  in  the  present  study.  At  each  adaptation 
step,  a  fixed  fraction  /  =  10%  of  the  candidate  elements  with  the  largest  error  indicators 
are  adapted.  The  marked  elements  are  refined  isotropically  through  its  local  mapping 
functions.  As  a  result,  the  newly  inserted  boundary  points  may  not  lie  on  the  geometry. 
In  order  to  ensure  the  accuracy  of  the  geometry  approximation,  all  of  the  newly  inserted 
boundary  points  are  remapped  to  the  real  geometry  by  querying  the  stored  geometry 
information  for  each  boundary  element.  Then  the  coordinates  of  the  interior  points  of  the 
modified  elements  are  updated  by  a  transfinite  interpolation  from  the  boundary  points. 

2.4  Numerical  Results 

2.4.1  Inviscid  Flow  over  the  NACA0012  Airfoil 

This  is  the  same  case  used  in  section  2.2.5.  To  assess  the  effectiveness  of  the  adaptation 
driven  by  the  adjoint-based  error  indicator,  h-adaptations  with  Cl  and  Cjj  as  the  output 
of  interest  are  performed  using  quadrilateral  elements  with  a  4th  order  scheme  (fc  =  3). 
The  newly  inserted  points  on  every  adaptation  stage  are  remapped  to  the  real  NACA0012 
airfoil  to  reduce  the  geometry  approximation  error.  The  initial  mesh  consists  of  560  p4 
curved  elements.  The  initial  mesh  and  the  final  adapted  meshes  using  lift  and  drag 
adjoints  are  shown  in  Figure  12.  Regions  near  the  trailing  edge  and  around  the  airfoil 
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(a)  The  initial  mesh 


(b)  Cl  Adjoint  (c)  C ;i  Adjoint 

Figure  12:  Adapted  mesh  for  the  adjoint-based  h-adaptation  for  a  inviscid  NACA  0012 
airfoil  at  Mq  =  0.5,  a  =  2°  (k  =  3) 


surface  are  refined  using  either  the  lift  or  drag  adjoint.  It  is  well-known  that  the  trailing 
edge  singularity  can  generate  spurious  entropy.  Therefore,  refinement  near  the  trailing 
edge  is  very  important  to  predict  an  accurate  drag  value.  Adaptations  using  lift-adjoint 
added  some  degrees  of  freedom  on  the  stagnation  streamlines,  where  the  lift  adjoint  is 
singular  and  oscillatory.  Figure  13  displays  the  Mach  and  adjount  solution  distributions  on 
the  initial  and  adapted  meshes.  Note  the  significant  improvement  in  solution  smoothness 
in  both  the  Mach  and  lift-adjoint  contours  from  the  initial  mesh  to  the  adapted  mesh. 
This  improvement  was  achieved  because  the  present  adaptation  framework  considers  both 
of  the  primal  and  adjoint  solutions. 

The  convergence  histories  of  the  lift  and  drag  coefficients  are  shown  in  Figure  14. 
The  corrected  outputs  are  computed  using  the  adjoint-based  error  estimate.  The  results 
show  that  the  corrected  values  converge  much  faster  than  the  uncorrected  ones,  and  all 
converge  to  the  same  value.  The  estimated  error  at  the  last  adaptation  stage  is  around 
10-10.  So  the  truth  Cl  and  Cd  hr  this  section  are  the  values  from  the  final  adapted  mesh. 
Figure  15  shows  the  output  error  of  the  adaptation  with  the  uniform  refinement  results 
for  comparison.  An  effective  convergence  rate  of  at  least  6th  order  was  achived  for  both 
Cl  and  Cd  with  h-adaptation.  It  is  clear  that  the  adjoint  based  h-adaptation  framework 
reduced  the  computational  cost  by  orders  of  magnitude  in  term  of  the  number  of  DOFs. 

2.4.2  Laminar  Flow  over  the  NACA0012  Airfoil 

In  this  case,  we  consider  subsonic  laminar  flow  over  a  NACA0012  airfoil  with  freestream 
Mo  =  0.5  and  angle  of  attack  a  =  1°.  The  Reynolds  number  based  on  the  chord  length 
of  the  airfoil  is  Be  =  5000.  The  same  initial  mesh  for  the  inviscid  test  case  was  used  in 
this  test.  The  drag  and  lift  coefficients  are  again  the  outputs  of  interest.  Adjoint-based 
h-adaptation  with  k  =  3  are  driven  by  the  output-based  error  indicator.  Additionally, 
uniform  h-refinement  is  performed  to  compare  those  adaptation  strategies. 

Figure  16  compares  the  Mach  contours  on  the  initial  mesh  and  the  final  adapted  mesh 


DISTRIBUTION  A:  Distributing  approved  for  public  release 


(a)  Mach  contours  on  the  initial  mesh 


(b)  Cl  adjoint  tj)2  on  the  initial  mesh 


(c)  Mach  contours  on  the  adapted  mesh 


(d)  Cl  adjoint  ip2  on  the  adapted  mesh 


Figure  13:  Primal  and  adjoint  solution  of  the  adjoint-based  h-adaptation  for  a  inviscid 
NACA  0012  airfoil  at  M0  =  0.5,  a  =  2°  (k  =  3) 


and  presents  the  adapted  meshes  from  different  adaptation  strategies.  Note  that  regions 
near  the  stagnation  streamlines  and  in  the  boundary  layer  were  targeted  for  refinements. 
The  trailing  edge  was  also  refined  repeatedly  to  reduce  the  effect  of  the  singularity. 

Figure  17  displays  the  convergence  of  the  lift  and  drag  coefficients  using  adjoint-based 
error  estimate  in  terms  of  nDOFs.  It  is  obvious  that  both  the  lift  and  drag  coefficients, 
from  all  adaptation  strategies,  converge  to  the  same  value.  Note  that  the  corrected 
outputs  converge  much  faster  than  the  uncorrected  ones.  The  truth  outputs  are  chosen 
from  the  output-based  h-adaptive  simulations  with  k  =  3,  whose  estimated  error  is  less 
than  10~8  at  the  final  stage.  Figure  18  compares  the  Cl  and  Cd  error  of  all  tested 
adaptation  strategies  with  results  from  uniform  h-refinements.  With  h-adaptation,  an 
effective  convergence  rate  of  6th  order  was  achieved  for  both  the  Cl  and  Cd ,  as  shown 
in  the  figure.  Again  it  is  shown  that  the  ajoint  based  h-adaptation  approach  can  reduce 
the  computational  cost  by  orders  of  magnitude. 

2.4.3  Inviscid  Flow  over  a  Sphere 

In  this  case,  we  consider  subsonic  inviscid  flow  over  a  sphere.  The  p3  hexahedral  mesh  is 
used  for  this  simulation.  Figure  19  shows  the  outline  of  the  computational  domain  and 
the  initial  surface  mesh  on  the  sphere.  The  inflow  Mach  number  is  set  to  be  0.3  with  an 
angle  of  attack  a  =  0°. 

First,  to  demonstrate  the  super  convergence  of  the  output  using  a  dual-consistent  CPR 
formulation,  a  uniform  refinement  study  is  performed.  As  shown  in  Figure  20,  a  super 
convergence  of  2k  + 1  is  obtained  for  the  Cd  error  with  k  =  1  and  k  =  2.  For  the  adaptive 
simulation,  a  relative  coarse  mesh  which  has  480  p3  hexahedral  elements  is  used  as  the 
initial  grid  and  the  3rd  order  CPR  scheme  with  the  Gauss  points  as  the  SPs/FPs  and  the 
LP  approach  is  tested.  The  adaptation  is  driven  by  the  adjoint-based  error  indicator  with 
drag  as  the  output  of  interest.  On  each  adaptation  level,  10%  of  the  current  elements  with 
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(a)  Cl  (b)  Cl  zoom  in 


Figure  14:  Cl  and  Cd  convergence  of  the  adjoint-based  h-adaptation  for  a  NACA  0012 
airfoil  at  Mq  =  0.5,  a  =  2°  (fc  =  3) 
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Figure  15:  Cl  and  Cd  error  of  the  adjoint-based  h-adaptation  for  a  NACA  0012  airfoil 
at  M0  =  0.5,  a  =  2°  (k  =  3) 


(a)  Mach  contours  on  the  initial  mesh  (b)  Mach  contours  on  the  adapted  mesh 


Figure  16:  Adjoint-based  h-adapted  mesh  for  a  NACA  0012  airfoil  at  Mq  =  0.5,  a  = 
1°, i?e  =  5000  (k  =  3) 
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Figure  17:  Cl  and  Cd  convergence  for  a  NACA  0012  airfoil  at  Mq  =  0.5,  a  =  1  °,Re  = 
5000  {k  =  3) 


(a)  Cl  Error 


(b)  Cd  Error 


Figure  18:  Cl  and  Cd  error  for  a  NACA  0012  airfoil  at  Mo  =  0.5,  a  =  1  °,Re  =  5000 
(k  =  3) 
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Figure  20:  Cd  error  for  the  inviscid  flow  over  a  sphere  problem. 


the  largest  error  are  marked  to  be  refined.  The  adapted  mesh  and  the  Mach  contours  are 
shown  in  Figure  21.  Regions  around  the  sphere  surface  are  refined  persistently.  Figure  22 
compares  the  drag  coefficient  error  of  the  adaptive  simulations  with  the  results  from  the 
uniform  h-refinements.  It  is  clear  to  see  that  the  current  adaptive  method  could  produce 
much  more  efficient  error  reductions  in  terms  of  the  number  of  the  DOFs.  An  effective 
convergence  order  of  6  is  obtained  through  the  adaptation,  which  is  much  faster  than  the 
uniform  refinements.  This  preliminary  adaptation  results  demonstrates  the  effectiveness 
of  the  present  adaptive  method  for  a  3D  problem. 


Figure  19:  The  initial  p3  curved  hexahedral  mesh  for  the  inviscid  flow  over  a  sphere 
problem. 
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(a)  Adaptation  level  1 


(b)  Adaptation  level  4 


Figure  21:  Mach  contours  on  the  adapted  mesh  for  the  inviscid  flow  over  a  sphere  at 
M0  =  0.3,  ct  =  2°  (k  =  2). 


Figure  22:  CD  errors  of  the  inviscid  flow  over  the  sphere  at  Mq  =  0.3,  a  =  2°(fc  =  2). 


2.4.4  Laminar  Flow  over  a  Sphere 

Next,  we  consider  steady  viscous  flow  over  a  sphere.  The  same  settings  from  the  Ref. 
[94,  95]  is  utilized  for  the  comparison  purpose.  The  Reynolds  number  based  on  the  sphere 
diameter  is  chosen  to  be  118.  The  inflow  Mach  number  is  0.2535  with  an  angle  of  attack 
a  =  0°.  The  3rd  order  CPR  scheme  with  the  Gauss  points  as  the  SPs/FPs  and  the 
LP  approach  is  tested.  Figure  23a  shows  the  initial  mesh  which  has  480  p3  hexahedral 
elements  and  the  corresponding  Mach  number  contours. 

The  adaptation  is  driven  by  the  adjoint-based  error  indicator  with  drag  as  the  output 
of  interest.  On  each  adaptation  level,  10%  of  the  current  elements  with  the  largest  error 
are  marked  to  be  refined.  The  adapted  mesh  and  the  Mach  contours  on  each  adaptation 
level  are  shown  in  Figure  23.  Regions  around  the  sphere  surface  are  refined  persistently. 
The  reference  Cjj  =  1.0162  is  chosen  from  Ref.  [95,  19].  Figure  24  compares  the  drag 
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coefficient  error  of  the  adaptation  with  the  result  from  the  uniform  h-refinements.  A 
super  convergence  rate  of  2k  is  obtained  for  the  uniform  h-refinement.  The  adjoint-based 
adaptation  shows  a  faster  convergence,  whose  effective  order  of  accuracy  is  5.9. 


(a)  The  initial  mesh  (b)  The  adapted  mesh,  level  3 

Figure  23:  Mach  contours  and  the  adapted  mesh  for  the  viscous  flow  over  a  sphere  problem 
at  M0  =  0.2535,  Re  =  118  and  a  =  0°  (fc  =  2). 


Figure  24:  CD  errors  of  the  viscous  flow  over  a  sphere  at  Mq  =  0.2535,  Re  =  118,  a  =  0° 
(k  =  2). 


2.4.5  Laminar  Flow  over  a  Delta  Wing 

In  this  test  case,  we  consider  laminar  flow  over  a  delta  wing  at  Mach  number  Mq  =  0.3 
with  a  high  angle  of  attack  a  =  12.5°.  This  case  is  a  part  of  the  International  workshop  on 
high-order  methods  (HOW)  and  the  EU  ADIGMA  project.  The  Reynolds  number  based 
on  the  root  chord  is  4000.  The  Prandtl  number  is  set  to  0.72  and  the  constant  viscosity 
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(a)  Adaptation  level  0  (b)  Adaptation  level  1  (c)  Adaptation  level  2 


(d)  Adaptation  level  3  (e)  Adaptation  level  4 


(f)  Adaptation  level  5 


Figure  25:  Mach  number  and  the  adapted  mesh  slices  on  different  adaptation  stages  for 
the  delta  wing  case. 


is  used.  The  delta  wing  has  a  swept  sharp  leading-edge  and  a  blunt  trailing  edge.  The 
Mach  number  contours  of  the  flow  around  the  delta  wing  at  different  adaptation  stages 
are  shown  in  Figure  25.  Both  of  the  singularity  along  the  leading  edge  and  the  regions 
around  the  smooth  vortices  are  targeted  to  refine.  The  Mach  number  distribution  is 
much  smoother  after  several  adaptation  levels.  Figure  26  shows  the  residual  history  of 
the  whole  adaptation  procedure.  For  each  adaptation  stage,  the  residual  norm  drops  10 
order.  Figure  27  compares  the  CD  error  and  the  CPU  time  of  the  adaptation  with  the 
result  from  the  uniform  refinements. 


2.4.6  Laminar  Flow  over  an  Analytic  3D  Body 

As  the  final  case,  we  consider  laminar  flow  over  a  streamlined  analytic  3D  body.  This  test 
case  is  a  part  of  the  International  workshop  on  high-order  methods  (HOW)  and  the  EU 
ADIGMA  project.  The  inflow  Mach  number  is  set  to  0.5  at  an  angle  of  attack  a  =  1°. 
The  Reynolds  number  is  5000  with  adiabatic  no-slip  wall  boundary  condition  enforced  on 
the  body  surface.  The  viscosity  is  assumed  to  be  a  constant  and  Prandtl  number  is  set 
to  0.72.  The  reference  area  is  0.05  and  the  reference  drag  coefficient  Cd  =  0.06317  from 
the  HOW  results  is  used.  The  coarsest  mesh  from  the  HOW  website  is  used  as  the  initial 
mesh,  which  consists  of  768  p4  curved  hexahedral  elements.  The  drag  adjoint  is  used  to 
drive  the  mesh  adaptation  procedures.  The  adapted  mesh  and  the  Mach  contours  are 
shown  in  Figure  28b.  The  refinements  are  mainly  performed  around  the  leading  edge  and 
around  the  body  surface.  Figure  29  display  the  residual  history  and  the  Cd  error. 
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Figure  26:  The  residual  history  of  the  whole  adaptation  procedure  for  the  delta  wing 
case. 


Figure  27:  Cd  error  and  the  CPU  time  of  the  delta  wing  case. 
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Res 


(a)  The  initial  mesh  and  Mach  contours 


(b)  The  adapted  mesh  and  Mach  contours 

Figure  28:  Laminar  flow  over  an  analytic  3D  body  at  M0  =  0.5,  a  =  1°  and  Re  =  5000. 


Figure  29:  Results  of  the  laminar  flow  over  an  analytic  3D  body. 
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Results  and  Discussion  on  Subgrid-scale  Stress  Model 
Evaluation 


3.1  Governing  Equation  and  SGS  Models 


The  governing  equations  for  three  dimensional  turbulent  flows  are  the  three  dimensional 
Navier-Stokes  equations.  As  a  simpler  counterpart,  we  consider  the  one  dimensional 
Burgers’  equation, 


du 

at 


du  d2u 

u^r  =  v^~ix  £ 
dx  ozx 


[-1.1], 


(36) 


where  u  is  the  state  variable  such  as  velocity,  v  is  a  constant  viscosity.  In  the  present 
study,  v  =  8 E  —  05  is  chosen  to  imitate  a  high  Reynolds  number  flow  problem.  To  derive 
the  LES  governing  equation,  we  apply  a  low-pass  spatial  filter,  G a{x,£)  satisfying  the 
following  conservative  property 


/OO 

GA(x,£)dZ=  1, 

-OO 


where  A  denotes  the  filter  width.  A  typical  filter,  the  box  filter,  is  defined  below, 


-£l< 


A 


=  2’ 
I  0  otherwise 


(37) 


(38) 


The  filtering  process  is  defined  mathematically  in  the  physical  space  as  a  convolution 
product.  The  filtered  variable  <j>(x,  t)  of  a  space-time  variable  <f>(x,  t)  in  ID  is  defined  as 


4>{x,t)=  /  G{x,£)<j)(x,t)dt. 


(39) 


The  filtering  process  is  linear,  i.e.  <j)  +  ip  =  </>  +  phi.  If  the  filter  width  is  constant,  the 
differential  and  the  filter  operators  commute,  i.e.  ^  In  the  present  study,  all 

the  filtering  processes  are  done  with  the  box  filter.  After  applying  the  filter  to  Eq.36,  we 


obtain 


du 

dt 


„ du  d2u 
U dx  d2x 


d{\uu  —  \ uu ) 
dx 


(40) 


The  unclosed  term  arises  due  to  the  filtering  of  the  nonlinear  convection  term 


tSGS 


-uu  —  -uu. 
2  2 


(41) 


This  is  the  SGS  of  the  Burgers’  equation.  SGS  models  act  as  the  closure  of  the  governing 
equation.  In  this  section  we  review  some  of  the  ideas  and  translate  them  to  work  for  the 
one  dimensional  Burgers’  equation. 


3.1.1  Static  Smagorinsky  Model 

The  SS  is  in  the  eddy  viscosity  form.  For  3D  incompressible  flow,  the  SGS  stress  is  defined 
as, 

TijGS  =  -ZvsGsSij,  (42) 

where  Stj  is  the  resolved  rate  of  strain  tensor,  and 

Si,j  =  \{dtUj  +djUi).  (43) 
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The  SGS  viscosity,  i'sgs  is  modeled  following  the  mixing  length  idea 


VSGS  =  (csA)2^2|^2,  (44) 

where  IS)2  =  SijSji,  cs  is  the  prescribed  coefficient.  By  comparing  the  mean  SGS  dissipa¬ 
tion  from  DNS  data  and  the  modeled  SGS  dissipation,  cs  can  be  determined.  Lilly  used 
this  procedure  for  isotropic  turbulence  to  obtain  cs  =  0.16.  The  SS  was  described  by  Moin 
and  Kim  Moin  P  [81],  Rogallo  and  Moin  Rogallo  RS  [86],  Lesieur  and  Metais  Lesieur  M 
[66]  and  Pope  SB  [88].  The  deficiency  of  this  model  first  showed  up  in  the  comparison 
of  the  modeled  SGS  stress  and  the  true  SGS  stress  computed  from  the  DNS  solution 
by  Clark  et  al.  Clark  RA  [22],  McMillan  and  Ferziger  McMillan  OJ  [79],  and  Bardina 
et  al.  Bardina  J  [6] .  The  comparisons  imply  that  the  model  does  not  capture  the  SGS 
adequately.  In  Meneveau  C  [80],  Meneveau  et  al.  gave  an  explanation  of  this  problem. 
Another  weakness  of  this  model  is  that  it  gives  non-zero  eddy  viscosity  in  laminar-flow 
regions.  Therefore  a  wall  function  is  needed  to  damp  the  SGS  viscosity  in  a  wall-bounded 
flow.  Next,  we  derive  its  ID  formulation.  The  rate  of  strain  in  ID  is 

"SGS  =  {csA)2\dxu\.  (45) 

Therefore  the  SGS  stress  becomes 

tsgs  =  _ UsGSs  (46) 

3.1.2  Dynamic  Smagorinsky  Model 

The  coefficient,  cs,  in  SS  is  prescribed.  However,  it  is  found  empirically  that  cs  depends 
on  the  flow,  being  0.1  for  plane  channel  flow  and  0.2  for  isotropic  turbulence  Durbin  P.  A 
[28].  The  DS  makes  it  a  variable  spatially  and  temporally.  It  introduces  a  test  filter  to 
the  resolved  scales  and  uses  the  assumption  of  scale  invariance  to  compute  the  model 
coefficient.  As  the  model  for  three  dimensional  turbulence  is  readily  available,  we  derive 
it  for  the  ID  Burgers’  equation  next. 

Following  Eq.  40,  we  consider  the  2nd  filter  with  width  A  ,  defined  as  A  =  qA.  By 
applying  this  filter  to  the  SGS  stress,  we  obtain 

j.SGS  _  _  i^nt.  (47) 

2  2 


By  applying  the  filter  to  the  LES  solution,  we  obtain  the  resolved  stress, 


r  1  1  -  - 
L  =  - mu  —  -uu. 

2  2 

(48) 

The  Germano  identity  can  be  written  as 

T  =  tsgs  +  L 

(49) 

where  T  =  \ uu  —  ^uu.  We  apply  the  SS  to  both  T  and  r  nd  assume 
coefficient,  cs, 

they  share  the  same 

-(csA)2 \dxu \dxu  =  ~(csA)2\dxu\dxw  +  L. 

(50) 

We  define 

M  =  A2 \dxu \dxm  —  A  2\dxu\dxu. 

(51) 

Thus  c2  =  jj.  It  is  assumed  that  cs  is  spatially  uniform  so  that  it  can  be  extracted  from 
the  test-filtering  operation  (Ghosal  et  al  1995)Ghosal  S  [39].  In  the  ID  test,  we  take  the 
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most  common  choice  of  7  =  2.  In  three  dimensions,  this  is  an  over-determined  system. 
To  minimize  the  square  error,  Lilly  used  the  following  approach 


2  _  (LjjMjj) 
s  {MijMnY 


(52) 


where  (•)  means  averaging  along  the  homogeneous  direction.  The  DS  gives  a  highly 
variable  eddy  viscosity  field  Germano  M  and  WH  [38]  including  negative  values  which 
makes  the  simulation  unstable.  Averaging  over  homogeneous  directions  was  used  by 
Germano  et  al.  Germano  M  and  WH  [38]  to  prevent  this  problem.  Ghosal  et  al.  Ghosal  S 
[39]  showed  that  this  procedure  minimizes  the  total  error  in  the  homogeneous  region  over 
which  the  averaging  is  performed.  With  these  modifications,  the  eddy  viscosity  still  can 
be  negative.  So  the  value  of  c2  is  clipped  to  be  non-negative.  In  ID,  we  don’t  have  these 
problems.  Thus  we  don’t  use  a  least  square  averaging  operation.  But  we  still  require  c2 
to  be  non-negative. 


3.1.3  Scale- Similarity  model 

The  SSM  was  first  introduced  by  Bardina  et  al.  Bardina  J  [6] .  It  assumes  scale  invariance 
between  the  computable  stress  L  and  the  SGS  stress  tsgs  .  This  assumption  was  verified 
with  empirical  band  pass-filtered  PIV  measurements  by  Liu  et  al.  Liu  S  [73].  It  suggests 
that  is  similar  to  a  stress  constructed  from  the  resolved  scales, 

tsgs  =  cssmL ,  (53) 

where  L  s  the  resolved  stress,  which  is  given  in  Eq.  48.  Many  different  second  filter  widths 
were  suggested  by  various  researchers.  The  Bardina’s  original  model  uses  the  same  filter 
width  for  the  two  filters,  i.e.  A  =  A  and  7=1.  Liu  et  al  used  7  =  2  and  Akhavan  et 
al  use  7=|  Akhavan  R  [2] .  The  coefficient  cssm  is  empirical  and  found  to  be  close  to 
1.  In  the  ID  test,  cssm  is  adjusted  to  be  0.25  with  ,  based  on  an  analysis  performed  in 
Z.J.  Wang  [120].  In  Bardina  J  [6],  the  true  and  modeled  stresses  showed  a  high  degree  of 
correlation  in  Bardina  et  al’s  a  priori  tests,  and  the  SSM  allowed  for  energy  backscatter. 
However,  this  model  was  found  to  be  not  sufficiently  dissipative.  Energy  accumulated  at 
small  scales  and  finally  led  to  numerical  instability.  In  the  present  study,  we  will  duplicate 
this  result  with  a  non-dissipative  numerical  scheme,  and  will  show  that  the  phenomenon 
does  not  occur  with  the  dissipative  FR/CPR  method. 


3.1.4  Mixed  Model 

To  resolve  the  above-mentioned  problem  of  the  SSM,  the  DS  is  included  in  the  formulation 
to  add  extra  dissipation.  In  three  dimensions,  the  mixed  model  (MM)  is 

TfjGS  =  CssmL ij  -  2 VSGsSij.  (54) 

Liu  et  al.  showed  that  the  magnitude  of  the  similarity  term  is  much  larger  than  that  of 
the  dissipative  DS  term.  Hence,  the  high  correlation  of  the  SSM  is  not  degraded  by  the 
extra  viscosity.  Zang  et  al.  used  this  model  for  recirculating  flows  with  7  =  1.  Wu  and 
Squires  applied  this  model  successfully  with  Lagrangian  averaging  in  simulations  of  3D 
boundary  layers  Wu  X  [114].  There  are  dynamic  ways  to  determine  cssm  as  well.  Vreman 
et  al  proposed  a  two-parameter  dynamic  MM  in  which  cs  and  cssm  are  both  calculated 
dynamically  with  7  =  1.  In  the  present  one  dimensional  study,  the  values  cssm  =  0.25 
and  7  =  2  are  used.  The  SGS  stress  is  defined  as 

T  =  Cssm  ( ~ UZt  ^uu)  -  VSGsS-  (55) 
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As  will  be  shown  later,  the  numerical  instability  of  the  SSM  model  is  not  a  problem  for 
numerical  schemes  with  embedded  numerical  dissipation.  The  study  will  be  described 
and  discussed  in  Section  4. 


3.1.5  Linear  Unified  RANS-LES  Model 


The  wall-bounded  turbulent  flows  at  high  Reynolds  number  are  a  significant  challenge  for 
LES.  The  near  wall  region  requires  a  high  resolution  grid  to  resolve  the  small  energetic 
scales.  The  linear  unified  RANS-LES  model  (LUM)  combines  RANS  with  LES  to  solve 
this  problem.  The  model  equations  for  incompressible  flows  are 

(  DUi  ,  nd(v+vt)Sik 


p  3  )  I  <~) 

dxi  dxk 


Dkt  _  +  Vtg2  _  o(l-co)fct 


^=Cul‘ZnutS2~%fl  u,2  + 


tl 

d(^+~n 


- ) _ Sm _ 1 

i  '  partialx  _■  / 


where  Cu i,  Cu 2,  Cu,  C/,  Grander u  are  all  model  constants,  Ut  is  the  filtered  velocity,  kt 
is  the  turbulent  kinetic  energy,  w  is  the  specific  dissipation,  77  is  the  time  scale  and  ut 
is  the  modeled  viscosity.  In  this  work,  we  only  focus  on  the  LES  aspect  of  the  model. 
Therefore  77  is  calculated  with  77,  =  lastA/k 2,  where  Z*  =  |.  In  summary,  for  the 
one-dimensional  Burgers’  equation,  we  solve 

{Du  _  9(^+^t)gg 
Dt  dx  gk 

Dkt  _  d((v+vt )-p)  ,  r,  /flu-12  _  o  (1-cq  )kt 

nt  dx  t^dx'  tl 


where  vt  =  tl  =  l^A/k^ 


3.2  Numeric  Methods 

3.2.1  The  High-order  CPR/FR  Method  for  the  ID  Burger’s  Equation 

To  give  a  complete  picture  of  the  LES,  in  this  section,  we  briefly  review  the  FR/CPR 
method’s  formula  for  the  ID  Burgers’  equation.  Huynh  Huynh  [54]  developed  a  high- 
order  FR/CPR  formulation,  which  was  later  employed  for  the  Navier-Stokes  equations 
on  hybrid  3D  meshes  T.  Haga  and  Wang  [94],  It  has  been  used  for  ID,  2D  and  3D 
laminar  and  turbulent  flows.  Validations  and  successful  applications  can  be  found  in  Gao 
and  Wang  [35]  Gao  and  Wang  [36]  Shi  and  Wang  [92]  Shi  and  Wang  [91]  Yu  et  al.  [118] 
Li  and  Wang  [69].  In  this  study,  we  apply  the  3rd  order  FR/CPR  scheme  to  the  ID 
Burgers’  equation  and  evaluate  its  performance  with  various  SGS  models.  The  FR/CPR 
formulation  for  the  inviscid  ID  Burgers’  equation  is  given  as 

+  A',.  I +  1/2 alpha^n  -  1/2)  =  0,  (56) 

where  utj  is  the  solution  at  solution  point  j  of  element  i,  Ui  is  the  solution  polynomial  for 
element  i,  n(  )  denotes  the  projected  flux  derivative  at  the  solution  point,  DeltaXi 
is  the  length  of  element  i,  [E"]i+ 1/2  and  [Fn]i_ i/2  are  the  differences  between  the  local 
flux  and  the  common  Riemann  flux  at  the  right  and  left  interfaces  of  element  i,  cxrj  and 
otL,j  are  the  correction  coefficients  independent  of  the  solution  variables.  For  the  viscous 
term  on  the  right  hand  side  of  Burgers’  equation,  we  follow  the  BR2  approach  Bassi  and 
Rebay  [10].  The  ID  version  is  described  below.  First  we  introduce  a  new  variable  R  = 

The  corrected  gradient  is  then 


Ri,j  — 


1  ^  /W  r  com  1  1  r  com  „  i  \ 

+  -^p(aRdlU  ~  ui\i+l/2  +  0!L,j[U  —  WiJi-1/2)) 
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where  [ucom]i+ 1/2  and  [ucom]i_1/2  are  the  common  solutions  at  interfaces,  [f(,:]i+i/2  and 
are  the  solutions  from  element  i  at  the  interfaces.  The  common  solution  is 
defined  as 


u 


com  _ 

*+1/2  — 


Ui+ 1/2 


+  U 


+ 

i+ 1/2 


2 


(58) 


where  ui+1j2  =  [u2:]j+i/2  and  uf+1,2  =  [wj+i]»+i/2  are  the  solutions  at  the  left  and  right 
sides  of  interface  i  +  1/2.  Next,  the  viscous  flux,  Fv  =  v^,  at  solution  points  can  be 
calculated  by 

-  /--"i /•’,,)•  (59) 

Then  can  be  obtained  by  using  the  Lagrange  polynomial  approach.  The  common 
viscous  flux  at  the  interface  is  needed  to  correct  at  solution  points, 


c)v 

rpv.com  _  i  ,(  \ com 

ri+ 1/2  “  ^QxJi+1/2- 


(60) 


For  the  common  gradient, 

(£)*+”  2  =  +  +  (£&i/2  (61) 

where  (^)~+1/2  and  (ur ^x])t+i/i  are  the  gradients  of  the  solution  of  the  left  and  right 
cells  with  no  correction,  r~+1j2  and  r^+1^2  are  the  corrections  to  the  gradients  due  to  the 
common  solution  at  the  interface.  More  specifically,  the  corrections  are, 


r2+l/2  =  7 7Z  ( Q.+  [uCOm 

where  alpha ~  and  a+  are  the  interface  correction  coefficients. 


]*+ 1/2)) 

(62) 

+]i+l/2), 

(63) 

3.2.2  Temporal  Discretization 


The  explicit  SSP  three-stage  3rd  order  Runge-Kutta  scheme  Gottlieb  and  Shu  [43]  is  used 
as  the  temporal  discretization.  Here  we  give  a  brief  description.  Rewrite  the  discretized 
Burgers’  equation  as 

dU 

—  =Res(U),Ut0  =  U0,  (64) 


where  Res{U)  is  a  function  of  solution  U  and  t.  Given  solution  Un,  we  obtain  solution 
Un+1  using 


[/«  =  un  +  AtRes(Un), 

(65) 

£/(  2)  =  hjn  +  Jf7(1)  +  tRes(U{1)), 

(66) 

Un+1  =  1Un  +  -t/(2)  +  2  A  tRes(U{2)). 

3  3  3 

(67) 

DISTRIBUTION  A:  Distributing  approved  for  public  release 


Figure  30:  Initial  condition  (left)  and  the  initial  energy  spectrum  (right) 


3.3  Initial  and  Boundary  Conditions  for  ID  Burgers’  Equation 


To  imitate  turbulence,  the  initial  energy  spectrum  is  given  in  the  Fourier  space  k.  In  the 
present  work,  the  following  initial  spectrum  is  used 


E0(k) 


A5~l  1  <  fc  <  5, 
Ak~i  k  >  5, 


(68) 


where  k  is  an  integer  varying  from  1  to  1280.  For  each  k,  the  velocity  u  has  a  random 
phase  angle,  /3,  in  [— 7r,7r]. 


n 

u(x)  =  y~^(2  E0(kj))i  sin(kjX  +  Pi)  +  1  (69) 

i 


A  is  a  constant  to  make  the  turbulence  intensity  A-  =  0.7%,  where  v!  =  -\A^'=1  W — ~ )> 
u  =  1 .  The  two  boundaries  are  set  to  be  periodic  due  to  the  periodicity  of  the  initial 
condition. 


3.4  Grid  and  Spatial  filter 

The  computational  domain  is  [—1, 1].  A  mesh  refinement  study  indicated  that  8,192  cells 
with  the  3rd  order  FR/CPR  method  are  required  to  resolve  all  the  scales.  Figure  31  shows 
the  energy  spectrum  at  for  the  linear  wave  propagation  with  the  same  initial  condition. 
There  is  no  visible  decay  at  even  the  highest  frequency.  We  consider  the  simulation  at 
this  resolution  a  DNS,  and  denote  A xdns  =  sl§2-  the  a  Pri°ri  study,  the  DNS  solution 
is  filtered  with  a  box  filter  with  A  =  32Axdns-  The  filtered  solution  is  obtained 

N  K+ 1 

'U’iyj  —  \  y  ^ 

n=  1  m—1 

where  N  is  the  number  of  cells  in  the  filtering  stencil  of  the  current  degree  of  freedom 
and  K  is  the  degree  of  the  polynomial  of  the  solution.  In  each  cell,  a  Gauss  quadrature 
rule  was  implemented  and  w  is  the  weighting  coefficient.  Then  the  filtered  solution  on  the 
DNS  grid  is  projected  to  the  (coarse)  LES  grid  if  necessary  to  serve  as  the  LES  solution. 
In  the  current  LES  cell,  the  projected  solution  at  each  solution  point  j  is  calculated  using 

r  K+ 1  N  K+l 

I  y  LjUi}jdx  —  y  (  /  Ljz  y  '  lj'U'7i,jdX)  (^1) 

J  Axles  j  n—1  ^  ^xdns  j 
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Figure  31:  Turbulent  energy  spectrum  at  t  =  0  and  t  =  T 
Table  5:  Correlation  coefficients  of  a  priori  test 


/\XLES 

Axons 

SS 

DS 

SSM 

Mix 

LUM 

1 

-0.10 

0.6 

0.95 

0.89 

-0.09 

2 

-0.10 

0.6 

0.95 

0.89 

-0.09 

4 

-0.10 

0.6 

0.95 

0.89 

-0.09 

8 

-0.10 

0.6 

0.95 

0.89 

-0.09 

16 

-0.10 

0.6 

0.95 

0.89 

-0.09 

32 

-0.09 

0.59 

0.95 

0.88 

0.04 

where  Lj  is  the  shape  function  defined  based  on  the  solution  points  of  the  LES  cell,  lj  is 
the  shape  function  based  on  the  solution  points  of  the  DNS  cell.  In  the  a  posteriori  study, 
we  do  the  same  thing  to  the  DNS  initial  condition  to  generate  the  LES  initial  condition. 

In  Figure  32,  A  =  8A xdns,  A xles  =  4A xdns  are  used  to  demonstrate  the  filtering 
operation.  Different  cell  sizes  for  LES  were  tested  to  evaluate  the  influence  of  the  trun¬ 
cation  error  and  the  SGS  modeling  error.  We  call  the  filter  used  on  the  initial  condition 
the  first  filter  and  the  filter  used  in  deciding  the  coefficient  of  the  dynamic  model  or  com¬ 
puting  the  resolved  SGS  stress  the  test  filter.  The  test  filter  width  is  2  times  the  width 
of  the  first  filter,  which  makes  7  =  2. 

3.5  Numerical  Results  and  Discussions 

In  this  section,  the  results  for  the  a  priori  and  a  posteriori  tests  are  presented.  Due  to 
the  nonlinear  convection  term,  shock  waves  start  to  appear  after  a  certain  time.  Thus 
all  results  are  obtained  at  a  time  T  =  0.1  ,  when  the  solution  is  still  smooth.  Figure  33 
shows  the  energy  spectrum  at  t  =  0  and  t  =  T  of  the  DNS.  We  can  see  that  the  high 
frequencies  are  damped  out  by  the  physical  viscosity  while  the  lower  frequencies  remain. 

3.5.1  A  Priori  Tests 

Figure  34  shows  the  SGS  stress  computed  using  different  models  based  on  the  filtered-DNS 
data  at  t  =  T  with  various  mesh  resolutions  and  a  fixed  filter  width  of  A  =  32  A  xdns- 
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■a -  DNS  solution 


Figure  32:  Comparison  of  various  solutions  (square:  fine  mesh  solution  points;  circle 
coarse  mesh  solution  points) 


Figure  33:  The  energy  spectrum  at  two  different  times 
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SGS  Stress  SGS  Stress  SGS  Stress 


True  SGS  stress 


Figure  34:  The  SGS  stress  comparison  in  the  a  priori  tests 
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Table  6:  Correlation  coefficients  of  a  posteriori  test 


Axles 
Ax  DNS 

SS 

DS 

SSM 

Mix 

LUM 

i 

-0.08 

0.6 

0.95 

0.89 

-0.06 

2 

-0.08 

0.6 

0.95 

0.89 

-0.07 

4 

-0.08 

0.6 

0.95 

0.89 

-0.07 

8 

-0.08 

0.6 

0.95 

0.89 

-0.07 

16 

-0.08 

0.6 

0.95 

0.89 

- 

32 

-0.10 

0.57 

0.92 

0.85 

- 

The  ratio  between  the  cell  size  of  LES  and  DNS  is  (a)  1,  (b)  2,  (c)  4,  (d)  8,  (e)  16,  (f)  32. 
For  the  SS  model,  cs  is  set  to  the  default  value  of  0.2  for  all  of  the  comparisons.  For  ILES, 
the  SGS  stress  is  0  everywhere.  From  Figure  34,  we  can  make  some  general  observation 
regardless  of  mesh  resolutions: 

•  No  models  are  able  to  predict  the  true  stress  in  both  amplitude  and  phase  (peaks 
and  valleys). 

•  Both  the  SSM  and  MM  always  correctly  predict  the  phase  of  the  true  stress. 

•  SS  correctly  predicts  the  phase  of  the  true  stress  about  half  the  time,  and  DS  agrees 
with  the  SS  when  the  phase  is  correct.  When  the  stress  computed  with  SS  has  a 
wrong  sign,  DS  sets  the  stress  to  0. 

•  LUM  agrees  very  well  with  SS  in  SGS  prediction  with  enough  grid  resolution,  but 
diverges  for  the  coarsest  mesh. 

Obviously  the  good  phase  prediction  capability  of  the  MM  is  due  to  the  dominant  SSM 
term.  Next  we  examine  the  correlation  of  the  modeled  stress  with  the  true  stress.  Table 
5  presents  the  correlation  coefficients  between  the  true  SGS  and  the  ones  computed  with 
the  models.  Clearly  the  SSM  and  the  MM  models  perform  the  best.  The  mesh  resolution 
Axles  does  not  have  any  significant  influence  on  the  model  behavior,  except  for  LUM. 
To  further  evaluate  the  behavior  of  these  models  in  an  actual  computation,  we  perform  a 
posteriori  tests  next. 

3.5.2  A  Posteriori  Tests 

In  this  test,  the  filtered  ID  Burgers’  equation  is  solved  with  different  models  on  different 
meshes  with  a  fixed  filter  size  A  =  32Axdns  ■  The  results  at  the  same  physical  time 
t  =  T  are  compared.  Figure  35  shows  the  SGS  stress  computed  using  different  models 
with  various  mesh  resolution. 

In  Figure  35,  we  can  see  that  the  results  are  very  similar  to  those  in  the  a  priori  test. 
We  can  draw  the  same  conclusions  here.  Table  6  shows  the  correlation  coefficients  for  all 
the  a  posteriori  tests.  The  SGS  stresses  computed  by  SSM  and  the  MM  always  show  high 
correlations  with  the  true  SGS  stress.  The  DS  comes  the  second.  The  SS  and  the  LUM 
models  yield  very  low  correlation  with  the  true  SGS  stress.  LUM  diverged  for  some  cases 
and  the  correlation  is  not  available.  The  LES  mesh  resolution  does  not  have  a  significant 
influence  on  the  model  behavior,  except  for  LUM. 

Figure  36  shows  the  comparison  of  the  solution,  u,  for  .  The  solution  computed 

with  the  true  SGS  stress  is  right  on  top  of  the  filtered  DNS  solution.  The  solutions 
computed  with  models  and  ILES  all  show  differences  with  the  filtered  DNS  solution. 
Table  7  shows  the  L2  norm  error  comparison.  When  the  LES  mesh  is  sufficiently  fine, 
it  is  clear  that  the  SSM  and  the  MM  produced  the  best  solutions.  This  is  because  both 
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SGS  Stress  SGS  Stress  SGS  Stress 


True  SGS  stress 


Figure  35:  The  SGS  stress  comparison  in  the  a  posteriori  tests 
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- B -  filtered  DNS  solution 

-  -Q - True  SGS  stress 


Figure  36:  Solution  Comparison 


Table  7:  L2  norm  error  of  the  solution 


&XLES 

Axdns 

True  stress 

ILES 

SS 

DS 

SSM 

Mix 

LUM 

1 

2.03E-08 

1.14E-05 

1.57E-05 

8.88E-06 

7.02E-06 

6.19E-06 

1.80E-05 

2 

2.13E-08 

1.14E-05 

1.57E-05 

8.88E-06 

7.02E-06 

6.18E-06 

3.02E-05 

4 

4.60E-08 

1.14E-05 

1.57E-05 

8.88E-06 

7.01E-06 

6.17E-06 

5.62E-05 

8 

4.91E-07 

1.14E-05 

1.57E-05 

8.97E-06 

6.94E-06 

6.25E-06 

1.07E-04 

16 

1.09E-05 

1.46E-05 

2.01E-05 

1.49E-05 

1.18E-05 

1.35E-05 

- 

32 

1.38E-04 

1.38E-04 

1.41E-04 

1.39E-04 

1.37E-04 

1.39E-04 

- 
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models  show  the  best  correlation  with  the  true  SGS  stress.  When  the  LES  mesh  is  coarse, 
the  truncation  error  is  dominant.  The  results  with  any  model  and  with  the  true  SGS  stress 
are  comparable.  Clearly  ILES  is  the  best  choice  because  it  costs  the  least. 

3.5.3  Sensitivity  of  the  Models  to  the  Mesh  Resolution 

Given  the  fixed  filter  width,  we  compare  models’  behavior  on  different  mesh  resolution. 
Figure  37  shows  different  modeled  SGS  comparison  with  respect  to  different.  The  ratio  = 
A xdns  '  can  see  ^at  both  in  the  a  priori  and  the  a  posteriori  tests,  all  the  models 

shows  no  sensitivity  to  Axles  except  for  the  LUM. 

3.5.4  Effects  of  Truncation  Error  vs.  SGS  Model  Error 

n  large  eddy  simulations,  the  numerical  results  depend  on  many  factors,  including  the 
flow  condition,  the  initial  and  boundary  conditions,  the  numerical  method,  the  computa¬ 
tional  mesh,  the  filter  and  the  SGS  model.  Some  of  the  factors  are  physical  and  others 
are  numerical,  and  they  intertwine  together  to  produce  the  final  solution.  At  the  most 
fundamental  level,  the  filter  width  A  in  a  LES  is  perhaps  the  most  critical  parameter, 
and  Pope  discussed  the  importance  of  the  filter  width  in  SB  [88] .  The  true  LES  solution 
can  be  obtained  by  filtering  the  DNS  solution  using  this  A.  In  reality,  however,  the  filter 
width  is  often  implicitly  tied  with  the  mesh  size.  In  such  cases,  mesh  refinement  conver¬ 
gence  studies  become  impossible  to  perform  because  the  filter  size  is  always  a  variable. 
One  can  only  see  convergence  when  the  mesh  size  approaches  that  required  of  a  DNS 
simulation.  Generally  speaking,  we  want  to  accurately  predict  the  SGS  stress  using  the 
numerical  solution  at  the  “resolved  scale”.  The  filtered  solution  u  is  always  taken  to  be 
the  solution  at  the  “resolved  scale”.  Let’s  consider  the  box  filter  here.  When  a  solution 
is  filtered  with  a  width  A  ,  we  often  state  that  waves  with  shorter  wavelengths  than  A 
are  filtered  out.  In  fact,  waves  of  wavelengths  of  2A  and  4A  are  heavily  damped  out  too. 
Based  on  our  analysis,  we  can  see  that  the  amplitudes  of  2A  and  4A  waves  are  reduced 
by  36%  and  10%  respectively  Z.J.  Wang  [120].  If  we  accept  36%  filtering  error  as  accept¬ 
able,  the  “resolved  scale”  should  be  2A  instead  of  A.  In  addition,  numerical  methods  also 
have  limited  resolution  depending  on  the  “points  per  wave”  (PPW)  or  “degrees  of  freedom 
per  wave”  (DOFPW).  Let’s  assume  that  for  the  present  3rd  order  FR/CPR  scheme,  9 
DOFPW  is  required  to  resolve  a  wave.  In  other  words,  3  elements  are  needed  for  a  wave 
since  there  are  3  DOFs  in  one  element.  A  truly  resolved  scale  must  meet  the  accuracy 
requirement  from  both  the  filtering  operator  and  the  numerical  scheme.  In  this  particular 
case,  the  resolved  scale  is 

SR  =  max(3 Axles,  2A),  (72) 

In  order  to  have  the  resolved  scale  determined  by  the  given  filter  width,  Axles  should 
satisfy  the  following  requirement 

2A 

A xLES  <  -j-  (73) 

In  the  case  of  second-order  finite  volume  methods,  each  element  has  1  solution  unknown. 
If  one  requires  20  PPW  for  accuracy,  the  resolved  scale  is  then 

Sr  =  max(2Q Axles,  2A).  (74) 

If  one  chooses  Axles  as  the  filter  width,  the  resolved  scale  is  20  times  larger  than  the 
filter  width  because  of  the  accuracy  requirement.  In  other  words,  the  numerical  truncation 
error  is  dominant  in  the  LES  results.  This  is  the  reason  why  we  see  smaller  and  smaller 
differences  between  the  ILES  and  LES  with  SGS  models  with  the  Axles  increase. 
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Figure  37:  The  modeled  SGS  with  different  mesh  resolution 
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4  Investigation  of  Scale  Similarity 

From  the  comparison  in  the  previous  section,  we  see  that  the  stress  computed  with  the 
scale  similarity  model  shows  the  highest  correlation  with  the  true  SGS  stress.  And  we 
discovered  that  they  are  related  by  a  factor  of  4,  i.e, 

—  =  4-  (75) 

1~true 

When  we  first  saw  the  relation  75,  we  suspected  that  this  was  a  bug.  After  we  tested 
different  initial  conditions,  we  were  convinced  that  it  was  not  a  bug.  The  fact  that  relation 
75  is  true  for  arbitrarily-generated  random  initial  conditions  prompted  us  to  look  for  a 
deeper  reason  resulting  in  the  following  analysis. 

4.1  Analysis  of  Scale  Similarity  with  a  Single  Fourier  Mode 

For  the  sake  of  simplicity  without  loss  of  generality,  we  consider  periodic  data  u(x)  at  a 
given  time  on  domain  [— 7r,  7t]  .  The  solution  can  be  decomposed  into  the  following  Fourier 
modes 

inf 

u{x)  =  ^2<inelnx,  (76) 

n—0 

where  i  =  \J—l,  and  n  is  the  wave  number.  To  illustrate  the  basic  idea,  we  first  consider 
a  single  Fourier  mode,  i.e.,  u(x)  =  emx  and  the  top  hat  filter.  The  filtered  solution  is 
then 

1  fx+ t  .  ,  nA 

Kx)  =  x  /  e”*dt  =  sinc(-)  ■  emx,  (77) 

a./*-  a  2 

where  sinc(^)  =  sa"A2  .  Obviously  the  filter  only  changes  the  magnitude  of  the  solution, 
but  not  the  phase.  In  addition,  we  have 

1  fx+ t 

uu( x)  =  -  /  ei2nid£  =  sincin A)  •  ei2nx.  (78) 

Ay*-# 

The  SGS  stress  is  then 

7,  *  *  •  /  a  \  i2nx  •  2  /  \  i2nx 

r  =  uu  —  uu  =  sinc{nA)  •  e  —  smc  (— — )  •  e 

T)  /\ 

=  [. sinc(nA )  —  smc2(— )]e'2ni.  (79) 

Next  we  apply  a  second  filter  with  a  width  of  A2  =  7A  to  the  resolved  variable  to 
obtain 

~  1  rx~\  2~  tlA  A 

u(x)  =  —  /  u(Odti  =  sinc(—)sinc(^—)emx,  (80) 

7A  yx-^  2  2 

and 

~  1  2  nA  ^/n A 

uifx)  =  -r  /  u{S)u{m  =  sinc2(  —  )sincC—)ei2nx.  (81) 

7A  Jx-i£  2  2 

The  SGS  stress  of  the  resolved  scale  is  then 

L  =  uu-  uu  =  sine2 (r^-)(sinc('ynA)  —  sinc2(^  9  ))el2nx .  (82) 
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From  Eqn.79  and  Eqn.  82,  we  obtain 


L  sine2  (^)  (sinciffnA)  —  sinc2{ -^^)) 

t  sinc(nA )  —  sinc2(^) 


In  the  limit  of  small  nA,  we  have 


sincinA)  =  1  —  +  0(nA)4. 

6 


Therefore,  we  obtain 


L  [1  -  ^  +  O(nA)4]  [-  +  0(?rA)4] 

-^  +  ^+0(nA)4 


=  72  +  O(nA)2.  (85) 


Note  that  the  error  term  is  quadratic.  In  the  special  case  of  7  =  2,  L  =  4 r.  As  it  turns 
out  this  result  is  also  true  for  the  Gaussian  filter.  The  filtered  solution  with  a  Gaussian 
filter  is 

/inf  6  0(a._£)2  0  /'inf  6(x—g)2 

\fi — — ^-)e  “a5  ein£dl ;  =  sqrt( — — =-)  /  e  ^  e?'n£d£.  (86) 

-inf  TrA  J  —  inf 

Set  A  =  ,  so  that  ^  =  x  — %A,  =  —A=dX.  Thus,  we  have 


n  /  6  A  /'IU1  — X2+m(a: — 4l)  ,y 

u(x)  =  \  /  e  T  V  ^  'dX 

V  ttA2  V6  7_i„f 

/"T"  -  /■“/  y2  .  A  y-  .  (tiA2 

=  J  —  etnx  e~x  ~znvsxdX  =  emxe~  — . 


V  PI  J-inf 

Similarly  we  can  derive  the  following  result 


„ . .  rni  /  6  6(,-o 

uu(x)=  I 


ei2na:d£ 


./-in/  V  n  aa- 

^  „i2nx  f  "~X‘ ^i2n~%X dX  =  ei2nx  ^ ^ 


\  PiC  J-inf  6 


The  SGS  stress  is  them 


T  =  uu-uu  =  ei2nxe - -  -  ei2nxe — =  ei2nxe - -  [1  -  e — n~].  (89 

Again  we  apply  a  second  filter  with  a  width  of  A2  =  7A  to  the  resolved  field  to  obtain 


//{./•)  =  e  24  //!./•)  =  e  ’  e 


(nA)2  (l.A)2 


V*  .  _  ("A)2  _  —  ( n  A)2  _(TnA)2  ■„ 

uyx)  —  e  12  uu(x)  =  e  12  e  6  e  . 


The  SGS  stress  of  the  resolved  scale  is 


(l  +  272)(-A)2  (-ynA)2 


L  =  iiu- uii  =  elInx  e  12  [1  —  e  12  ]. 
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From  89  and  92,  we  obtain 


(7nA)2 

L  (l-272)(nA)2  1  —  e  12 


=  e 


(n  A)2 

1  —  e  12 


In  the  limit  of  small  nA,  we  have 


(93) 


—  =  (1  +  O(nA)2) 1  _  1  _  +  O(nA)4 - 2 


1  -  1  - 


(nA)2 

12 


0(nA)4 


=  Y  +  0(nA)2. 


(94) 


4.2  Analysis  of  Scale  Similarity  with  All  Fourier  Modes 

Next  we  consider  a  solution  with  all  the  Fourier  modes,  i.e. , 

inf 

u(x)  =  55  anemx . 

n— 0 

With  the  top  hat  filter,  we  obtain  the  following  filtered  solution 


(95) 


inf 


nA , 


i(x)  =  55  ansinc(—)  ■  emx. 


n— 0 


In  addition,  we  have 


uu(x)  =  £  f]  anamMC(n  +  m)Aei("+’")l 


n— 0  m—0 


The  SGS  stress  is  then 


(96) 


(97) 


inf  inf 


-  -  -  r  •  (n  +  ra)A  .  .nA  .  .raA., 

=  uu  —  uu  =  2_^  cLnam[sinc - - - sinc{——  )sinc (— — )] e  v  ^  ;  . 


n=0  m— 0 

Next  we  apply  a  second  filter  with  a  width  of  7A  to  the  resolved  variable 

inf 

n.„  sinrt  - 

2  ’  v  2 


\  '  ■  /  tiA  .  .  .  yu  A  ^nx 

=  2  yanStnc{  —  )sinc(—^—)e  , 


n=0 


(98) 


(99) 


and 


V'  V'  •  +  j(n+r; 

=  2^  anamsmc(  —  )smc(——)stnc( - - - )e  v 


um 


n—0  m— 0 

The  SGS  stress  of  the  resolved  scale  is  then 


(100) 


L  =  uu  —  wu 

st'  st'  ■  ,nA  .  mA  .  'y(n  +  m)A 

=  2^  2^  anamSinc(  —  )smc(^-)(smc( - - - ) 

n—0  m—0 

-  sincC^)sinc(^^))e^n+m)x . 

Now  let’s  consider  each  term  in  Eqn.  98  and  101.  It  is  obvious  that 


(101) 
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Lnm  sinc( ^ )sinc{  )(sinc(  T(ra+m) A )  _  sinc(  ) sinc(  'y™A ) ) 

Tnm  sinc^n+^^A  —  sinc{r!YL)sinc(Jn£L) 

= 'y2  +  0[(n  +  m)A]2 .  (102) 

Therefore,  we  have 

—  =  72  +  0[{n  +  m)A]2.  (103) 

T 

n  the  same  limit.  The  analysis  with  the  Gaussian  filter  is  similar  and  is  not  repeated 
here. 


4.3  Analysis  of  Scale  Similarity  in  2D 

In  two  dimensions,  we  only  perform  a  single  mode  analysis  with  the  top  hat  filter.  Consider 
the  following  two  dimensional  velocity  field 

u(x,  y )  =  einxeimv,  v(x,  y)  =  eipxeiqy.  (104) 


The  filtered  solution  is  then 


r*+  §  rv+ § 


nA, 


,m  A, 


eir*eimr,dtndt  =  eim:eimy  •  sinc{  —  )sinc{— ), 


JV~~2 


'j{x,y)  = 


1 


A2 

In  addition,  we  have 


rx+ 1  /•"  •  2  .  ,  .  nA  oA 

I  /  eip^ ezqri dr/dt;  =  eipxeiqy  ■  sinc(—)sinc{  —  ). 

a-#  2  2 


(105) 

(106) 


l 

Mx,y)  =  ^  A 

^  X  2 

=  ei(n+p)xei(m+q)y  .  j  (  (»  +P)A  g  (»*  +  g)A 
v  2  2 

The  SGS  stress  is  then 


rv+  2 


’y-f 


^+P)iei(m+q)rldT}d^ 


(107) 


r  =  uv  —  uv 
=  ei(n+p)xei(m+q)y  [sinc( 


(n  +  p)  A  .  {m  +  q)  A 
- - - )sinc{ - - - ) 


/  uA  _  .  77? A .  .  .pA.  .  ,gA 

—  sz?rc(  -^-)sinc{  —^—)sinc{-^-)sinc(  )J . 

Applying  a  second  filter  with  a  width  7  A  to  the  resolved  variable,  we  obtain 

~  inx  imv  .  ,nA  .  ,mA  .  ,7 nA  .  7mA 

u  =  e  e  y  ■  sinc(——)sinc(——)sinc(  )sinc{ — - — ), 


and 


~  ivx  iav  ■  ■  ,7PA  ■  ,79^ 

v  =  ep  e  qv  ■  sinc(  —  )sinc(  —  )smc{ - )sinc( - ), 

1  2  2  2  v  2 


Denote  a  =  sinc(n^)sinc(ln^)sinc(^-)sinc(^).  Then  we  have 


2  2 

hi)  =  auv,  uv  =  auv. 


(108) 

(109) 

(HO) 

(HI) 
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The  SGS  stress  of  the  resolved  scale  is  then 


L  =  uv—  uv  =  a(uv  —  uv) 

=  r,An+^Am^v[sincC{n  +  p)A  WW7(m  +  g)  A 


)sinc(- 


2  '  y  2 

7n  A  .  ,7mA  .  ,7pA  .  ,7gA 


-) 


-)sinc(- 


2  >inc(— )si»c(— )]. 


From  108  and  112,  we  obtain 


L  sincC(n+oP)A)sincC(m+q)A)  -  sinc(^)sinc(^)sinc{^)sinc{^) 


—  =  a- 

T 


(112) 


(113) 


sinc(  ^n+^A)sinc(  (m+g)A )  _  Sjnc(i^)Sjnc(22^)sjnc(£A)Sjnc(2A) 

In  the  limit  of  small  (n  +  m  +  p  +  q)  A,  we  have  alpha  «  1,  and 

CD 

sinc(A  +  B)  ■  sinc{C  +  D)  —  sine  A  ■  sincB  ■  sincC  ■  sincD  = - - - 1-  HOT.  (114) 

O 

Finally,  we  derive  the  following  result  using  Eqn.  114 

L 


l2- 


(115) 


4.4  Implications  for  Large  Eddy  Simulation 

The  present  analysis  shows  that  perfect  scale  similarity  exists  for  arbitrary  (periodic) 
data  including  turbulence  under  the  assumption  that  the  spectrum  contains  relatively 
low  frequency  contents  with  respect  to  the  filter  width,  regardless  of  amplitude  and  phase 
angle  of  each  mode.  Obviously  for  an  arbitrary  spectrum  including  both  high  and  low 
frequency  contents,  the  present  analysis  is  not  valid.  This  is  easily  seen  in  Figure  38, 
which  displays  the  modeled  and  true  SGS  stress  based  on  the  full  spectrum  shown  in 
Figure  30,  using  the  same  filter  width  which  is  16 A dns-  The  correlation  between  the 
modeled  and  true  stresses  is  quite  low. 

Next  let’s  examine  whether  Eqn.  115  is  true  in  an  actual  LES.  The  promise  of  the 
SSM  is  that  the  SGS  stress  is  highly  correlated  with  the  stress  computed  based  on  the 
resolved  scale,  taken  to  be  u.  Take  the  top  hat  filter  for  example.  Modes  of  smaller 
wavelength  than  D  corresponding  to  the  cutoff  wavenumber  k/±  are  filtered  out.  In  LES, 
it  is  believed  that  the  SGS  stress  from  higher  modes  close  to  the  cutoff  wave  number  k\ 
plays  an  important  role.  In  the  next  test,  we  therefore  include  modes  between  and  2&a 
using  a  filter  width  D/2  to  filter  the  spectrum  shown  in  Figure  30.  The  filtered  solution  is 
then  treated  as  DNS  data,  which  is  used  to  obtain  the  true  stress.  This  true  stress  is  also 
compared  with  the  stresses  computed  using  the  SSM  based  on  the  resolved  scale,  i.e.,  u. 
Two  test  filter  widths  are  used  corresponding  to  7  =  1  and  2.  The  results  are  displayed  in 
Figure  38.  Note  that  there  is  a  reasonably  high  level  of  correlation  between  the  stresses. 

The  ratio  between  the  true  and  modeled  stresses  are  computed  using  simple  averages 


L=  (L) 
r  (r)  ' 


(116) 


The  correlation  coefficients  and  the  average  stress  ratios  from  10  realizations  are  summa¬ 
rized  in  Table  5.  The  table  confirms  that  the  true  stress  shows  a  quite  high  correlation 
with  the  modeled  stress,  with  an  average  correlation  coefficients  of  0.88  and  0.69  for  7  =  1 
and  2,  respectively.  In  addition,  7  =  1  demonstrates  consistently  higher  correlation  coef¬ 
ficients  than  7  =  1.  This  may  indicate  that  one  should  use  the  same  filter  width  for  the 
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Figure  38:  The  true  stress  and  the  modeled  stress  for  the  full  spectrum 


- B -  True  stress 

—  -A-  -  L  with  gamma=1 


Figure  39:  Comparison  between  the  true  stress  computed  with  the  SGS  between  and 
2 and  the  modeled  stresses  computed  using  u 
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second  filter  in  an  SSM  implementation.  Furthermore,  the  ratio  of  the  averaged  stresses 
remains  a  constant  with  different  realizations,  indicating  that  this  ratio  is  only  dependent 
on  the  spectrum.  However,  the  ratio  is  much  smaller  than  7.  This  result  appears  to 
agree  well  with  others  in  the  literature  AW  [4]  Liu  S  [73]. 

4.5  Investigation  of  Stability  of  Scale  Similarity  Model 

In  Bardina’s  original  paper,  the  SSM  was  found  unstable  in  some  simulations  when  used 
with  a  central  finite  difference  scheme.  To  remedy  the  instability,  a  MM  with  the  DS 
model  was  developed  to  stabilize  the  simulations.  In  this  section,  we  attempt  to  show 
that  the  extra  dissipation  added  by  the  MM  is  not  necessary  for  the  FR/CPR  method 
which  has  embedded  numerical  dissipation  to  automatically  damp  high  frequency  modes. 
We  first  demonstrate  that  there  is  indeed  a  pile-up  of  high  frequency  modes  with  a  central 
difference  scheme  in  solving  nonlinear  equations  such  as  the  Burgers’  equation,  while  there 
is  no  such  pile-up  with  a  dissipative  high-order  FR/CPR  scheme.  For  this  purpose,  we 
conduct  a  numerical  study  with  the  initial  condition  of  a  single  Fourier  mode, 

u(x)  =  2(£/(l))1  /2  sin(TTx)  +  1.  (117) 

where  =  2(£'o(l))1^2  =  0.012.  The  ID  inviscid  Burgers’  equation  is  employed  to  mimic 
very  high  Reynolds  number  problems.  We  run  the  simulation  until  t  =  26  when  it  is  right 
before  a  shock  wave  develops.  First,  the  upwind  flux  and  the  central  flux  are  employed  in 
the  3rd  order  FR/CPR  scheme  to  compare  their  behaviors.  Figure  40  shows  the  energy 
spectrum  at  t  =  26  with  different  mesh  resolutions.  On  the  finest  mesh,  both  the  central 
and  upwind  schemes  produced  a  converged  solution  within  the  visible  energy  spectrum  in 
the  figure.  On  the  two  coarser  meshes,  we  can  see  clearly  that  energy  is  piling  up  at  high 
frequencies  on  those  meshes  for  the  simulation  with  the  central  flux.  But  the  upwind 
flux  is  able  to  smoothly  damp  out  the  high  frequency  modes  so  that  they  are  never 
accumulated  to  cause  stability  problems.  Next  we  test  the  influence  of  the  SGS  models 
on  the  energy  spectrum.  Figure  41  shows  the  spectrum  comparison  of  the  simulations 
with  and  without  the  SSM  and  MM.  The  filter  width  equals  to  the  cell  size.  We  can  see 
that  with  the  central  flux,  the  SSM  neither  damps  out  all  the  energy  accumulated  at  high 
frequencies  nor  accumulates  more  energy  there.  Thus  the  extra  dissipation,  i.e.  the  DS, 
is  necessary  to  stabilize  the  simulation.  It  is  worth  noting  that  the  extra  dissipation,  in 
the  MM,  also  damps  out  the  energy  at  some  lower  frequencies,  which  does  harm  to  the 
resolved  large  scales. 

We  also  verify  that  a  central  difference  finite  difference  scheme  behaves  similarly  with 
the  CPR  scheme  with  a  central  flux.  Figure  42  indeed  shows  that  the  4th  order  central 
finite  difference  method  has  a  similar  performance  to  the  3rd  order  FR/CPR  scheme  with 
the  central  flux.  This  means  that  for  schemes  that  are  not  dissipative,  more  dissipation 
may  be  necessary  to  stabilize  the  turbulent  flow  simulations  with  the  SSM  model.  But 
for  the  dissipative  ones,  such  as  the  FR/CPR  method  with  an  upwind  flux,  no  extra 
dissipation  is  needed. 
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- B -  CPR-3  upwind  24  dots 

-  -O - CPR-3  central  24  dots 


k 


Figure  40:  The  spectrum  of  the  upwind  flux  and  the  central  flux  with  3rd  order  FR/CPR 
scheme  at  t  =  26 
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Figure  41:  The  spectrum  of  different  models  with  3rd 


order  FR/CPR  scheme  at  t  =  26 
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Figure  42:  The  spectrum  of  the  3rd  order  FR/CPR  scheme  and  the  Ath  order  finite 
difference  scheme  with  the  central  flux  at  t  =  26 


5  Results  and  Discussion  on  MeshCurve:  An  Auto¬ 
mated  Low- Order  to  High- Order  Mesh  Converter 

5.1  Low  versus  High-Order  Meshes:  An  Illustrative  Example 

Figure  43  illustrates  the  difference  between  a  low-order  mesh  and  a  high-order  mesh  with 
a  side-by-side  comparison  of  two  circular  meshes.  The  key  difference  is  the  number  of 
nodes,  with  the  low-order  mesh,  shown  on  the  left,  having  fewer  than  the  high-order 
mesh,  shown  on  the  right.  The  extra  nodes  enhance  simulation  accuracy  by  serving  as 
secondary  interpolation  points.  Also,  they  enhance  geometric  accuracy  by  tracing  the  arc 
of  curved  edges. 

The  goal  of  meshCurve  is  to  transform  a  mesh  similar  to  the  one  on  the  left  of  figure 
43  into  a  mesh  similar  to  the  one  on  the  right,  but  in  3D.  Figure  44  is  a  before-and-after 
image  from  meshCurve.  The  additional  nodes  are  not  visible  in  the  high-order  mesh,  but 
the  surface  curvature  is  clearly  discernible. 

5.2  The  Design  of  meshCurve 

5.2.1  Feature  Criteria 

From  our  requirements,  we  assembled  the  following  list  of  features  for  meshCurve. 

•  Full  support  for  3D  unstructured  CGNS  meshes — any  cell  type  in  any  combination. 

•  Support  for  multi-zone,  multi-patch  CGNS  meshes,  including  the  ability  to  selec¬ 
tively  reconstruct  interior  patches  without  affecting  far-field  boundaries. 

•  Mesh  processing  without  CAD  geometry  files. 

•  Mesh  processing  in  the  face  of  complex  edge  arrangements  and  sharp  corners. 


DISTRIBUTION  A:  Distributing  approved  for  public  release 


(a)  a  low-order  mesh  (b)  a  high-order  mesh 

Figure  43:  A  low-order  mesh  contrasted  with  a  high-order  mesh.  The  key  difference  is 
that  the  high-order  mesh  has  more  nodes  than  the  low-order  mesh.  The  extra  nodes 
enhance  the  geometric  accuracy  by  tracing  the  path  of  curved  edges. 


(a)  a  low-order  mesh  (b)  a  high-order  mesh 

Figure  44:  Before-and-after  comparison  of  a  mesh  upgraded  with  meshCurve.  The  addi¬ 
tional  nodes  in  the  high-order  mesh  are  not  visible  but  the  surface  curvature  is  visible. 
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•  Easy-to-use  graphical  interface,  requiring  minimal  effort  to  learn  and  navigate. 

•  Interactive  3D  graphics  to  show  the  mesh  and  to  provide  visual  feedback  after  the 
upgrade. 

•  Cross-platform  support — Linux,  Windows,  and  Mac. 

•  Solid  code  base  with  minimal  bugs,  designed  for  maintainability  and  future  up¬ 
grades. 

•  Reasonably  low  memory  footprint  and  fast  operation  on  a  desktop  computer. 

•  Minimal  reliance  on  outside  software  libraries  to  increase  control  of  the  code  and 
simplify  design. 

5.2.2  Design  Decisions 

To  satisfy  requirements,  we  made  the  following  decisions. 

•  Programming  would  be  in  CH — h  due  to  the  speed  and  flexibility  of  the  language. 
We  opted  for  C++11,  the  most  recent  version  of  the  language,  because  of  newly 
introduced  optimizations  and  platform-independent  multithreading.  With  multi¬ 
threading,  it  would  be  possible  to  utilize  the  multiple  cores  of  modern  CPUs. 

•  To  support  our  desired  cross-platform  graphical-user-interface  (GUI),  we  selected 
the  popular  Qt  application  framework.1  The  Qt  framework  is  a  platform-independent 
windowing  library  which  includes  multiple  user  interface  widgets  and  support  li¬ 
braries.  Applications  built  with  Qt  adapt  to  the  target  operating  system.  Conse¬ 
quently,  meshCurve  appears  to  be  an  XI 1  application  when  run  on  Linux,  a  Win¬ 
dows  application  when  run  on  Microsoft’s  OS,  and  a  Macintosh  application  when 
run  on  OS  X.  Figure  45  shows  an  annotated  view  of  the  meshCurve  user  interface 
as  it  appears  on  Microsoft  Windows.  The  figures  elsewhere  show  the  interface  as  it 
appears  on  Linux  (Ubuntu)  and  Mac  OS  X. 

•  For  interactive  3D  visualization,  we  chose  the  Visualization  Toolkit  (VTK),  an  open 
source,  cross-platform  3D  visualization  system  created  by  Kitware.[89]  VTK  wraps 
the  underlying  graphics  APIs,  simplifying  cross-platform  development  while  provid¬ 
ing  a  capable  visualization  environment. 

•  For  fast  linear  algebra,  we  decided  to  use  the  high-speed  Armadillo  C++  library.  [87] 
Armadillo  interfaces  platform-specific  LAPACK  and  BLAS  libraries  to  provide  vec¬ 
torized  mathematics  and  fast  matrix  routines.  Its  syntax  is  based  on  Matlab,  making 
mathematical  programming  both  simple  and  efficient. 

•  We  constructed  our  own  mesh  database.  Building  our  own  database  provided  free¬ 
dom  to  optimize  for  specific  needs,  while  also  granting  flexibility  to  design  for  future 
extensibility. 

5.3  User  Workflow 

Figure  46  diagrams  the  mesh-processing  workflow  from  the  user’s  perspective.  The  first 

step  of  the  workflow  is  the  loading  of  a  mesh  file,  either  CGNS  or  Gmsh  format.  The 

loaded  mesh  displays  as  a  3D  view-  rotatable,  zoomable,  and  translatable  via  the  mouse. 

The  view  can  optionally  be  set  to  wireframe  or  opaque  solid,  and  separate  patches  within 

xhttp:/ /www.qt.io 


DISTRIBUTION  A:  Distributing  approved  for  public  release 


Figure  45:  The  meshCurve  graphical  interface,  as  it  appears  on  Microsoft  Windows  7. 
By  default,  three  panels  are  visible:  1)  3D  viewing,  2)  configuration  setting,  and  3) 
information  display.  The  information  panel  and  the  configuration  panel  can  be  hidden  to 
increase  the  size  of  the  3D  viewport.  A  row  of  toolbar  buttons  and  menus  along  the  top 
provide  options  to  the  user. 


Figure  46:  The  mesh  upgrade  workflow  from  the  user’s  perspective. 
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the  mesh  can  be  drawn  individually.  Also,  surface  nodes  can  be  displayed  on  the  mesh. 
Viewing  the  nodes  is  especially  useful  as  a  final  verification  after  the  upgrade  has  com¬ 
pleted. 

The  next  step  is  to  select 
patches  to  process.  During  the 
upgrade,  the  selected  patches  will 
be  given  curvature.  For  efficiency, 
unselected  patches  will  not  be 
given  curvature  although  they 
will  be  given  high-order  nodes. 

Figure  47  illustrates  patch  selec¬ 
tion,  showing  a  fighter  jet  mesh 
with  wing  and  tail  patches  se¬ 
lected.  Green  highlight  marks  the 
selection. 

Next  comes  the  identification 
of  feature  curves,  such  as  the  rim 
of  an  airplane  canopy  or  the  sharp 
edge  of  a  box.  Feature  curves  de¬ 
fine  the  geometry,  so  it  is  impor¬ 
tant  that  meshCurve  know  the  lo¬ 
cation  of  the  curves  within  the 
mesh.  Every  feature  curve  may 
be  manually  flagged  edge-by-edge,  via  the  mouse.  Alternatively,  the  curves  can  be  flagged 
automatically  by  way  of  meshCurve ’s  own  feature-curve  detection  engine.  Employing 
discontinuity  detection  algorithms  developed  by  Jiao  and  Bayyana,[59]  the  feature-curve 
detection  engine  can  identify  curves  even  when  they  cross  nearly  flat  portions  of  the  mesh. 
After  auto-detection,  the  mouse  can  be  used  to  fine  tune  the  result. 

Once  feature  curves  have  been  flagged,  upgrading  can  commence.  During  the  upgrade, 
surface  geometry  is  approximated  by  a  continuous  2D  polynomial,  which  is  used  to  posi¬ 
tion  high-order  nodes  on  the  surface,  giving  the  surface  curvature.  If  the  curvature  is  too 
great,  the  cells  abutting  the  boundary  will  intersect.  To  prevent  intersections,  meshCurve 
can  bend  the  edges  of  interior  cells  to  coincide  with  the  boundary.  This  bending  action 
is  applied  by  default  but  may  be  turned  off  by  a  checkbox  in  the  configuration  panel  for 
surface  reconstruction,  as  a  way  to  speed  calculation.  Not  all  meshes  need  to  have  their 
interior  deformed  in  this  way. 

At  the  completion  of  the  upgrade,  the  3D  view  updates  to  show  the  new  high-order 
mesh.  The  user  can  then  save  the  mesh  as  a  new  CGNS  file. 

5.4  Case  Studies 

Three  case  studies  are  presented  to  illustrate  the  capability  of  meshCurve.  For  each  case 
study,  a  mesh  is  upgraded  from  linear  to  2nd  order  (i.e.  low  to  high-order).  A  picture 
of  the  mesh  is  presented,  rendered  with  meshCurve.  A  table  follows,  summarizing  key 
information  to  demonstrate  performance.  The  table  lists: 

1.  the  total  number  of  cells 

2.  the  total  number  of  nodes 

3.  the  fraction  of  patches  given  curvature 

4.  the  time  it  took  to  process  the  mesh 


Figure  47:  Screenshot  of  the  program  meshCurve,  il¬ 
lustrating  patch  selection.  The  green  portions  of  the 
mesh  are  the  patches  that  have  been  selected  for  pro¬ 
cessing. 
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5.  whether  or  not  interior  deformation  was  activated 


6.  minimum  and  maximum  Jacobians  of  the  original  low-order  mesh 

7.  minimum  and  maximum  Jacobians  of  the  new  high-order  mesh 

Processing  was  done  on  a  Macintosh  computer  with  a  2.7  GHz  Intel  i5  processor  (quad 
core). 

The  Jacobians  are  a  check  for  the  presence  of  intersecting  and/or  inverted  cells,  both 
of  which  could  be  problematic  for  simulations.  Negative  Jacobians  indicate  trouble.  Only 
one  of  the  case  studies  encountered  a  negative  Jacobian,  and  this  occurrence  was  reme¬ 
died  by  interior  deformation.  For  the  other  case  studies,  interior  deformation  had  little 
affect.  Since  the  timing  calculations  show  that  interior  deformation  can  be  the  most  time 
consuming  step,  it  is  best  to  apply  interior  deformation  only  when  it  is  truly  needed. 
meshCurve  applies  it  by  default,  assuming  the  worst,  but  interior  deformation  can  be 
easily  unselected.  It  should  be  noted  that  the  interior  deformation  algorithm  isn’t  guar¬ 
anteed  to  prevent  negative  Jacobians;  it  merely  reduces  their  likelihood.  For  example, 
negative  Jacobians  can  occur  if  boundary  cells  are  borderline  degenerate.  The  most  re¬ 
cent  version  of  meshCurve  has  the  ability  to  calculate  the  Jacobians,  alerting  the  user  to 
potential  problems. 
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Figure  48:  Case  study:  meshCurve  was  used  to  process  the  above  mesh,  upgrading  it  from 
linear  to  2nd  order.  The  mesh  represents  a  sphere  inside  a  square  shaped  domain.  There 
are  two  patches  in  the  mesh,  one  representing  the  surface  of  the  sphere  (highlighted  green) 
and  the  second  representing  the  outer  surface  of  the  domain  (shown  but  not  highlighted) . 
The  green  areas  of  the  mesh  were  given  curvature  during  the  upgrade.  Table  8  presents 
statistics  for  the  before  and  after  mesh. 


Table  8:  Statistics  for  figure  48  mesh:  meshCurve  was  used  to  upgrade  the  originally  linear  mesh  to  2nd 
order.  The  program  was  executed  on  a  Macintosh  computer  with  a  2.7  GHz  Intel  Core  i5  processor  running 


Mac  OS  X. 


#  cells 

#  nodes 

patches: 
recon. /total 

interior 

defer. 

time 

original  Jacobian 
[  low  —  high  ] 

new  Jacobian 
[  low  —  high  ] 

480 

588 

1/2 

no 

0.08  sec 

[  0.0131343  -  192.609  ] 

[  0.0129203  -  192.613  ] 

480 

588 

1/2 

yes 

0.11  sec 

[  0.0131343  -  192.609  ] 

[  0.0129203  -  192.613  j 
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5.5  Internal  Architecture  Overview 

The  mesh  database  is  the  core  of 
meshCurve.  It  manages  all  informa¬ 
tion  associated  with  the  mesh,  serving 
as  both  the  information  storage  facility 
and  the  conduit  through  which  the  var¬ 
ious  components  of  meshCurve  commu¬ 
nicate.  There  are  four  external  classes 
which  interface  the  database  to  provide 
additional  functionality.  These  classes 
perform  the  following  tasks:  1)  file  in¬ 
put/output,  2)  feature  curve  detection, 

3)  surface  reconstruction/interior  defor¬ 
mation,  and  4)  3D  visualization.  The 
code  for  the  user  interface  wraps  around 
the  database  and  accompanying  classes. 

It  exposes  the  internal  switches  for  the 
user  to  manipulate.  The  modular  struc¬ 
ture  of  the  code  base  and  its  object- 
oriented  design  make  the  program  easy 
to  edit.  Figure  51  illustrates  the  flow  of 
data. 

The  mesh  database  is  built  for  maximum  flexibility.  It  is  designed  with  no  assumption 
about  the  cell  types  or  node  orderings  stored  within.  These  definitions  are  placed  else¬ 
where,  as  a  configuration  table,  which  makes  updating  the  definitions  extremely  easy.  A 
second  example  of  flexibility  is  the  storage  scheme  for  element  connectivities,  which  uses 
both  forward  and  backward  references  to  not  only  answer  questions  such  as  “Which  nodes 
occupy  this  cell?”  but  also  the  inverse,  “Which  cells  contain  this  node?”  Furthermore,  the 
database  provides  a  mechanism  to  cycle  through  elements  as  a  list,  as  if  the  data  were 
stored  as  an  array.  The  database  also  provides  a  mechanism  to  jump  through  elements 
neighbor-to-neighbor,  as  if  the  mesh  were  stored  as  a  connected  graph. 

5.5.1  Key  Algorithms:  feature  curve  detection 

To  find  feature  curves,  meshCurve  employs  an  algorithm  proposed  by  Jiao  and  Bayyana, 
which  identifies  the  curves  probabilistically,  using  a  combination  of  several  topological 
metrics. [59]  The  algorithm  proceeds  through  several  steps.  Surface  nodes  and  edges  are 
flagged  as  candidate  discontinuities  if  the  local  dihedral  angle,  ridge  direction,  and  angle 
defect  fall  within  a  predefined  range.  Then  the  candidate  discontinuities  are  linked  into 
chains.  The  list  of  chains  is  filtered,  and  each  chain  receives  a  classification  based  on 
the  properties  of  its  starting  and  ending  nodes.  As  the  algorithm  progresses,  the  list  is 
narrowed  until  only  the  most  likely  candidates  remain.  In  the  end,  the  remaining  chains 
are  classified  as  feature  curves.  Figure  52  illustrates  the  capabilities  of  the  method, 
implemented  in  meshCurve.  Blue  lines  represent  detected  curves;  orange  dots  represent 
detected  corners. 

5.5.2  Key  Algorithms:  surface  reconstruction 

For  high-order  surface  reconstruction,  meshCurve  employs  a  variation  of  Weighted  Aver¬ 
aging  of  Local  Fittings  (WALF),  developed  by  Jiao  and  Wang.  [60]  WALF  reconstructs  a 
mesh  surface  by  averaging  a  set  of  locally  fitted  2D  Taylor  polynomials.  In  the  original 
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Figure  51:  The  internal  data  flow  for  the  pro¬ 
gram  meshCurve.  The  database  is  the  core 
of  the  program,  serving  as  storage  facility  and 
communication  hub.  Four  external  classes  in¬ 
terface  the  database  to  accomplish  the  four 
tasks  1)  File  input/output,  2)  3D  visualization, 
3)  feature  curve  detection,  and  4)  surface  recon¬ 
struction. 
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Figure  49:  Case  study:  meshCurve  was  used  to  process  the  above  mesh,  upgrading  it 
from  linear  to  2nd  order.  The  mesh  represents  a  static  mixer.  There  are  four  patches  in 
the  mesh,  one  of  which  is  shown  (highlighted  green).  The  other  three  patches  have  been 
hidden  from  view  to  make  the  visualization  easier  to  interpret.  Curvature  imparted  to 
the  high-order  mesh  is  clearly  discernible  around  the  rim  of  the  inlet  spout.  Table  ?? 
presents  statistics  for  the  before  and  after  mesh. 

Table  9:  Statistics  for  figure  49  mesh:  meshCurve  was  used  to  upgrade  the  originally  linear  mesh 
to  2nd  order.  The  program  was  executed  on  a  Macintosh  computer  with  a  2.7  GHz  Intel  Core  i5 
processor  running  Mac  OS  X. 


#  cells 

#  nodes 

patches: 

interior 

time 

original  Jacobian 

new  Jacobian 

recon. /total 

defor. 

[  low  —  high  ] 

[  low  —  high  ] 

13,761 

2,786 

4/4 

no 

15.3  sec 

[  0.00164224  -  0.0418735  ] 

[  0.000572914  -  0.043354  ] 

13,761 

2,786 

4/4 

yes 

16.1  sec 

[  0.00164224  -  0.0418735  ] 

[  0.000572914  -  0.0429421  ] 

Table  10:  Statistics  for  figure  50  mesh,  representing  an  airfoil:  meshCurve  was  used  to  upgrade  the 
originally  linear  mesh  to  2nd  order.  The  program  was  executed  on  a  Macintosh  computer  with  a 
2.7  GHz  Intel  Core  i5  processor  running  Mac  OS  X. 


#  cells 

#  nodes 

patches: 

interior 

time 

original  Jacobian 

new  Jacobian 

recon. /total 

defor. 

[  low  — 

high  ] 

[  low  — 

high  ] 

455,988 

478,840 

1  /  7 

no 

14.2  sec 

[  3.00E-016  - 

1.31E-008  ] 

[  -4.54E-018  - 

1.31E-008  ] 

455,988 

478,840 

1  /  7 

yes 

6.6  min 

[  3.00E-016  - 

1.31E-008  ] 

[  3.00E-016  - 

1.31E-008  ] 
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Figure  50:  Case  study:  meshCurve  was  used  to  process  the  above  mesh,  upgrading  it 
from  linear  to  2nd  order.  The  mesh  represents  an  airfoil.  There  are  seven  patches  in  the 
mesh,  one  representing  the  surface  of  the  airfoil  (highlighted  green)  and  six  representing 
the  outer  surfaces  of  the  simulation  domain.  In  the  above  view,  the  front  patch  has  been 
hidden  to  reveal  the  interior  airfoil.  Table  ??  presents  statistics  for  the  before  and  after 
mesh. 
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(a)  (b) 

Figure  52:  Screenshots  of  the  program  meshCurve,  illustrating  automatic  feature  curve 
detection.  Blue  lines  mark  feature  curves  identified  by  the  program.  In  the  left  subfigure, 
the  edges  of  a  spring  have  been  identified  as  feature  curves.  In  the  right  subfigure,  feature 
detection  has  been  limited  to  the  portion  of  the  mesh  isolated  for  processing.  The  mesh 
represents  a  fighter  jet  with  a  single  wing  selected.  Curve  detection  has  marked  the  wing’s 
boundary  as  a  feature  curve.  Interior  curves  and  corners  have  also  been  identified. 


algorithm,  the  polynomials  are  fit  to  every  node  on  the  surface.  In  our  version,  the  poly¬ 
nomials  are  fit  to  every  face.  The  use  of  faces  as  the  location-of-fit  reduces  the  required 
number  of  averages.  Also  the  use  of  face-centered-fits  facilitates  the  treatment  of  feature 
curves  as  boundaries,  thereby  preventing  the  feature  curves  from  being  smoothed  away. 
An  iterative  solver,  initialized  to  zero,  helps  reduce  oscillations. 

5.5.3  Key  Algorithms:  interior  deformation 

To  deform  the  interior  elements,  meshCurve  employs  an  interpolation  scheme  developed 
by  Luke,  Collins,  and  Blades  originally  designed  for  the  simulation  of  fluid  structure 
interaction. [76]  According  to  this  method,  each  interior  node  moves  in  accordance  with  a 
weighted  average  of  the  boundary  node  motion.  Every  boundary  node  is  weighted  inverse 
to  its  distance  and  in  proportion  to  the  area  it  commands. 
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6  Conclusions 


In  the  present  project,  a  robust  and  efficient  high-order  CFD  software  suite  for  the  com¬ 
pressible  Navier-Stokes  equations  that  can  provide  engineering  accuracy  for  real  world 
industry  problems  is  developed.  The  following  points  summarize  the  key  accomplish¬ 
ments  and  significant  conclusions  that  can  be  drawn  from  the  current  program: 

To  estimate  the  error  of  an  engineering  output,  we  extended  the  dual-weighted  residual 
method  originally  developed  in  the  variational  framework  to  the  high-order  CPR  method 
which  is  in  the  differential  form.  A  dual-consistent  CPR  formulation  of  hyperbolic  conser¬ 
vation  laws  is  developed  and  its  dual  consistency  is  analyzed.  Super-convergent  functional 
and  error  estimate  for  the  output  with  the  CPR  method  are  obtained.  Factors  affecting 
the  dual  consistency,  such  as  the  solution  point  distribution,  correction  functions,  bound¬ 
ary  conditions  and  the  discretization  approach  for  the  non-linear  flux  divergence  term, 
are  studied. 

Next,  we  developed  a  parallel  adjoint-based  adaptive  CPR  solver  with  the  capability  of 
handling  any  element-based  error  estimate  and  arbitrary  discretization  orders  for  mixed 
grids.  The  current  method  have  been  applied  to  aerodynamic  flows  and  challenging 
engineering  applications.  Numerical  tests  show  that  significant  savings  in  the  number  of 
DOFs  can  be  achieved  through  the  adjoint-based  adaptation. 

Five  SGS  models  are  evaluated  with  the  ID  Burgers’  equation  discretized  with  the 
CPR  method.  Different  LES  cell  sizes  were  tested  with  a  fixed  filter  width.  In  both 
the  a  priori  and  a  posteriori  tests  on  a  fine  LES  mesh,  the  SSM  and  the  MM  showed 
excellent  correlation  with  the  true  SGS,  while  the  other  models  do  not  predict  the  SGS 
stress  satisfactorily.  However,  as  the  LES  cell  size  increases,  numerical  truncation  error 
is  dominant  in  the  results.  In  this  case,  none  of  the  models  shows  any  benefits  over  ILES. 
The  analysis  of  scale  similarity  shows  that  perfect  scale  similarity  exists  for  arbitrary 
(periodic)  data  including  turbulence  under  the  assumption  that  the  spectrum  contains 
relatively  low  frequency  contents  with  respect  to  the  filter  width,  regardless  of  amplitude 
and  phase  angle  of  each  mode.  In  an  actual  large  eddy  simulation,  in  which  both  large 
and  sub-grid  scales  exist,  the  present  result  on  the  ratio  of  the  resolved  scale  stress  and  the 
SGS  stress  may  be  the  upper  limit.  Test  results  with  data  including  higher  modes  near  the 
grid  cutoff  demonstrate  that  there  is  a  high  level  of  correlation  between  the  modelled  and 
SGS  stresses.  Furthermore,  7  =  1  demonstrates  consistently  higher  correlation  coefficients 
than  7  =  2.  This  may  indicate  that  7  =  1  is  preferred  in  a  SSM  implementation.  The 
stability  of  the  SSM  is  also  investigated.  The  study  shows  that  it  is  the  central  flux  rather 
than  the  SSM  that  causes  energy  accumulating  at  high  frequencies,  which  may  lead  to 
the  instability  of  a  simulation.  In  this  case,  extra  dissipation  other  than  SSM,  the  DS,  for 
example,  is  necessary.  However,  the  schemes  with  upwind  flux  smoothly  damps  out  the 
energy  at  high  frequencies.  Thus  no  extra  dissipation  is  needed  to  stabilize  the  simulation 
with  the  SSM. 

Finally,  we  created  a  user-friendly  GUI-based  software  tool,  called  meshCurve,  to 
convert  linear  unstructured  meshes  to  curved  high-order  meshes.  Based  on  a  geometry 
reconstruction  procedure,  the  software  infers  boundary  curvature  from  the  mesh  itself, 
without  need  of  a  CAD  geometry  file.  The  software  aims  to  be  a  cross-platform  tool  for 
mesh  upgrading,  useful  wherever  high-order  meshes  are  needed. 
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With  increased  computational  power  and  progress  in  numerical  methods  over  the  past  several  decades, 
Computational  Fluid  Dynamics  (CFD)  is  now  used  routinely  as  a  powerful  tool  in  the  design  of  aircraft. 
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potential  has  been  demonstrated  conclusively  for  smooth  problems  in  the  latest  International  Workshops 
on  High-Order  Methods.  For  non-smooth  problems,  solution  based  hp-adaptation  offers  the  best  promise. 
The  primary  objective  of  the  present  study  is  to  develop  robust  and  efficient  high-order  CFD  methods  and 
tools  for  the  compressible  Navier-Stokes  equations  that  can  provide  engineering  accuracy  for  real  world 
industry  problems.  Several  pacing  items  are  addressed,  which  include  hp-adaptations,  sub-grid  stress 
models  for  large  eddy  simulations,  and  high-order  mesh  generation. 
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